transfer.build_matrices(mgdof);
MGLevelObject<Vector<double> > u(0, tr.n_levels()-1);
- MGTools::reinit_vector(mgdof, u);
+ reinit_vector(mgdof, u);
// First prolongate the constant
// vector. For Lagrange elements,
// the values are just the number
// $Id$
// Version: $Name$
//
-// Copyright (C) 2000 - 2007 by the deal.II authors
+// Copyright (C) 2000 - 2007, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <iomanip>
#include <iomanip>
#include <algorithm>
+#include <numeric>
using namespace std;
+template <int dim, typename number, int spacedim>
+void
+reinit_vector_by_blocks (
+ const dealii::MGDoFHandler<dim,spacedim> &mg_dof,
+ MGLevelObject<BlockVector<number> > &v,
+ const std::vector<bool> &sel,
+ std::vector<std::vector<unsigned int> >& ndofs)
+{
+ std::vector<bool> selected=sel;
+ // Compute the number of blocks needed
+ const unsigned int n_selected
+ = std::accumulate(selected.begin(),
+ selected.end(),
+ 0U);
+
+ if (ndofs.size() == 0)
+ {
+ std::vector<std::vector<unsigned int> >
+ new_dofs(mg_dof.get_tria().n_levels(),
+ std::vector<unsigned int>(selected.size()));
+ std::swap(ndofs, new_dofs);
+ MGTools::count_dofs_per_block (mg_dof, ndofs);
+ }
+
+ for (unsigned int level=v.get_minlevel();
+ level<=v.get_maxlevel();++level)
+ {
+ v[level].reinit(n_selected, 0);
+ unsigned int k=0;
+ for (unsigned int i=0;i<selected.size() && (k<v[level].n_blocks());++i)
+ {
+ if (selected[i])
+ {
+ v[level].block(k++).reinit(ndofs[level][i]);
+ }
+ v[level].collect_sizes();
+ }
+ }
+}
+
template <int dim>
void check_block(const FiniteElement<dim>& fe,
const vector<bool>& selected,
if (selected[i])
deallog << ' ' << i;
deallog << std::endl;
-
+
Triangulation<dim> tr;
GridGenerator::hyper_cube(tr);
tr.refine_global(2);
-
+
MGDoFHandler<dim> mgdof(tr);
DoFHandler<dim>& dof=mgdof;
mgdof.distribute_dofs(fe);
DoFRenumbering::component_wise(mgdof);
vector<unsigned int> ndofs(fe.n_blocks());
DoFTools::count_dofs_per_block(mgdof, ndofs);
-
+
for (unsigned int l=0;l<tr.n_levels();++l)
DoFRenumbering::component_wise(mgdof, l);
- std::vector<std::vector<unsigned int> > mg_ndofs(mgdof.get_tria().n_levels());
+ std::vector<std::vector<unsigned int> > mg_ndofs(mgdof.get_tria().n_levels(),
+ std::vector<unsigned int>(fe.n_blocks()));
MGTools::count_dofs_per_block(mgdof, mg_ndofs);
deallog << "Global dofs:";
for (unsigned int i=0;i<mg_ndofs[l].size();++i)
deallog << ' ' << mg_ndofs[l][i];
deallog << std::endl;
- }
-
+ }
+
PrimitiveVectorMemory<Vector<double> > mem;
MGTransferBlock<double> transfer;
transfer.build_matrices(dof, mgdof, selected);
transfer.initialize(factors, mem);
MGLevelObject< BlockVector<double> > u(0, tr.n_levels()-1);
-
- MGTools::reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
+
+ reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
// Prolongate a constant function
// twice
for (unsigned int b=0;b<u[2].n_blocks();++b)
deallog << '\t' << (int) (u[2].block(b)*u[2].block(b)+.4);
deallog << std::endl;
-
+
u[1] = 0.;
u[0] = 0.;
transfer.restrict_and_add(2,u[1],u[2]);
for (unsigned int b=0;b<u[0].n_blocks();++b)
deallog << '\t' << (int) (u[0].block(b)*u[0].block(b)+.5);
deallog << std::endl;
-
+
// Check copy to mg and back
// Fill a global vector by counting
// from one up
v.reinit (ndofs);
for (unsigned int i=0;i<v.size();++i)
v(i) = i+1;
-
+
// See what part gets copied to mg
u.resize(0, tr.n_levels()-1);
- MGTools::reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
-
+ reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
+
transfer.copy_to_mg(mgdof, u, v);
for (unsigned int i=0; i<u[2].size();++i)
deallog << ' ' << (int) u[2](i);
deallog.attach(logfile);
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
-
+
std::vector<double> factors;
FE_DGQ<2> q0(0);
FE_DGQ<2> q1(1);
FE_RaviartThomasNodal<2> rt0(0);
FE_RaviartThomasNodal<2> rt1(1);
-
+
FESystem<2> fe0(rt1, 1, q1, 1);
FESystem<2> fe1(rt0, 2, q0, 2);
-
+
vector<bool> s1(2, true);
deallog << "All" << std::endl;
check_block(fe0, s1, factors);
-
+
s1[1] = false;
deallog << "Velocity" << std::endl;
check_block(fe0, s1, factors);
-
+
s1[1] = true;
s1[0] = false;
deallog << "Pressure" << std::endl;
// $Id$
// Version: $Name$
//
-// Copyright (C) 2000 - 2007 by the deal.II authors
+// Copyright (C) 2000 - 2007, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <iomanip>
#include <iomanip>
#include <algorithm>
+#include <numeric>
using namespace std;
+template <int dim, typename number, int spacedim>
+void
+reinit_vector_by_blocks (
+ const dealii::MGDoFHandler<dim,spacedim> &mg_dof,
+ MGLevelObject<BlockVector<number> > &v,
+ const std::vector<bool> &sel,
+ std::vector<std::vector<unsigned int> >& ndofs)
+{
+ std::vector<bool> selected=sel;
+ // Compute the number of blocks needed
+ const unsigned int n_selected
+ = std::accumulate(selected.begin(),
+ selected.end(),
+ 0U);
+
+ if (ndofs.size() == 0)
+ {
+ std::vector<std::vector<unsigned int> >
+ new_dofs(mg_dof.get_tria().n_levels(),
+ std::vector<unsigned int>(selected.size()));
+ std::swap(ndofs, new_dofs);
+ MGTools::count_dofs_per_block (mg_dof, ndofs);
+ }
+
+ for (unsigned int level=v.get_minlevel();
+ level<=v.get_maxlevel();++level)
+ {
+ v[level].reinit(n_selected, 0);
+ unsigned int k=0;
+ for (unsigned int i=0;i<selected.size() && (k<v[level].n_blocks());++i)
+ {
+ if (selected[i])
+ {
+ v[level].block(k++).reinit(ndofs[level][i]);
+ }
+ v[level].collect_sizes();
+ }
+ }
+}
+
+
template <int dim>
void check_select(const FiniteElement<dim>& fe, unsigned int selected)
{
deallog << fe.get_name()
<< " select " << selected << std::endl;
-
+
Triangulation<dim> tr;
GridGenerator::hyper_cube(tr);
tr.refine_global(2);
-
+
MGDoFHandler<dim> mgdof(tr);
DoFHandler<dim>& dof=mgdof;
mgdof.distribute_dofs(fe);
DoFRenumbering::component_wise(mgdof);
vector<unsigned int> ndofs(fe.n_blocks());
DoFTools::count_dofs_per_block(mgdof, ndofs);
-
+
for (unsigned int l=0;l<tr.n_levels();++l)
DoFRenumbering::component_wise(mgdof, l);
std::vector<std::vector<unsigned int> > mg_ndofs(mgdof.get_tria().n_levels());
for (unsigned int i=0;i<mg_ndofs[l].size();++i)
deallog << ' ' << mg_ndofs[l][i];
deallog << std::endl;
- }
-
+ }
+
MGTransferBlockSelect<double> transfer;
transfer.build_matrices(dof, mgdof, selected);
// constant one and the l2-norm is
// the number of degrees of freedom.
MGLevelObject< Vector<double> > u(0, tr.n_levels()-1);
-
- MGTools::reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
+
+ reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
u[0] = 1;
transfer.prolongate(1,u[1],u[0]);
transfer.restrict_and_add(1,u[0],u[1]);
deallog << "u1\t" << (int) (u[1]*u[1]+.5) << std::endl
<< "u0\t" << (int) (u[0]*u[0]+.5) << std::endl;
-
+
// Check copy to mg and back
// Fill a global vector by counting
// from one up
v.reinit (ndofs);
for (unsigned int i=0;i<v.size();++i)
v(i) = i+1;
-
+
// See what part gets copied to mg
u.resize(0, tr.n_levels()-1);
- MGTools::reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
-
+ reinit_vector_by_blocks(mgdof, u, selected, mg_ndofs);
+
transfer.copy_to_mg(mgdof, u, v);
for (unsigned int i=0; i<u[2].size();++i)
deallog << ' ' << (int) u[2](i);
deallog.attach(logfile);
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
-
+
FE_DGQ<2> q0(0);
FE_DGQ<2> q1(1);
FE_RaviartThomasNodal<2> rt0(0);
FE_RaviartThomasNodal<2> rt1(1);
-
+
FESystem<2> fe0(rt1, 1, q1, 1);
FESystem<2> fe1(rt0, 2, q0, 2);
-
+
check_select(fe0, 0);
check_select(fe0, 1);
-
+
check_select(fe1, 0);
check_select(fe1, 1);
check_select(fe1, 2);
// $Id$
// Version: $Name$
//
-// Copyright (C) 2000 - 2007 by the deal.II authors
+// Copyright (C) 2000 - 2007, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <iomanip>
#include <iomanip>
#include <algorithm>
+#include <numeric>
using namespace std;
+
+template <int dim, typename number, int spacedim>
+void
+reinit_vector_by_blocks (
+ const dealii::MGDoFHandler<dim,spacedim> &mg_dof,
+ MGLevelObject<BlockVector<number> > &v,
+ const std::vector<bool> &sel,
+ std::vector<std::vector<unsigned int> >& ndofs)
+{
+ std::vector<bool> selected=sel;
+ // Compute the number of blocks needed
+ const unsigned int n_selected
+ = std::accumulate(selected.begin(),
+ selected.end(),
+ 0U);
+
+ if (ndofs.size() == 0)
+ {
+ std::vector<std::vector<unsigned int> >
+ new_dofs(mg_dof.get_tria().n_levels(),
+ std::vector<unsigned int>(selected.size()));
+ std::swap(ndofs, new_dofs);
+ MGTools::count_dofs_per_block (mg_dof, ndofs);
+ }
+
+ for (unsigned int level=v.get_minlevel();
+ level<=v.get_maxlevel();++level)
+ {
+ v[level].reinit(n_selected, 0);
+ unsigned int k=0;
+ for (unsigned int i=0;i<selected.size() && (k<v[level].n_blocks());++i)
+ {
+ if (selected[i])
+ {
+ v[level].block(k++).reinit(ndofs[level][i]);
+ }
+ v[level].collect_sizes();
+ }
+ }
+}
+
+
template <int dim>
void check_block(const FiniteElement<dim>& fe)
{
deallog << fe.get_name() << std::endl;
std::vector<bool> selected(fe.n_blocks(), true);
-
+
Triangulation<dim> tr;
GridGenerator::hyper_cube(tr);
tr.refine_global(2);
-
+
MGDoFHandler<dim> mgdof(tr);
DoFHandler<dim>& dof=mgdof;
mgdof.distribute_dofs(fe);
-
+
// Make sure all orderings are the
// same
Point<dim> direction;
DoFTools::count_dofs_per_block(mgdof, ndofs);
std::vector<std::vector<unsigned int> > mg_ndofs(mgdof.get_tria().n_levels());
MGTools::count_dofs_per_block(mgdof, mg_ndofs);
-
+
MGTransferPrebuilt<BlockVector<double> > transfer;
MGTransferBlock<double> transfer_block;
MGTransferBlockSelect<double> transfer_select;
transfer.build_matrices(mgdof);
transfer_block.build_matrices(dof, mgdof, selected);
transfer_select.build_matrices(dof, mgdof, 0);
-
+
BlockVector<double> u2(mg_ndofs[2]);
BlockVector<double> u1(mg_ndofs[1]);
BlockVector<double> u0(mg_ndofs[0]);
<< ' ' << v1.block(0).l2_norm()
<< ' ' << v2.block(0).l2_norm()
<< std::endl;
-
+
v2 = u2;
u1 = 0.;
u0 = 0.;
deallog << "Select " << v0.block(0).l2_norm()
<< ' ' << v1.block(0).l2_norm()
<< std::endl;
-
-
+
+
// Check copy to mg and back
BlockVector<double> u;
u.reinit (ndofs);
v.resize(0, tr.n_levels()-1);
wb.resize(0, tr.n_levels()-1);
ws.resize(0, tr.n_levels()-1);
- MGTools::reinit_vector_by_blocks(mgdof, v, selected, cached_sizes);
- MGTools::reinit_vector_by_blocks(mgdof, wb, selected, cached_sizes);
- MGTools::reinit_vector_by_blocks(mgdof, ws, 0, cached_sizes);
+ reinit_vector_by_blocks(mgdof, v, selected, cached_sizes);
+ reinit_vector_by_blocks(mgdof, wb, selected, cached_sizes);
+ reinit_vector_by_blocks(mgdof, ws, 0, cached_sizes);
deallog << "copy to mg";
transfer.copy_to_mg(mgdof, v, u);
transfer_block.copy_to_mg(mgdof, wb, u);
transfer_select.copy_to_mg(mgdof, ws, u);
-
+
for (unsigned int l=0; l<tr.n_levels();++l)
{
wb[l].add(-1., v[l]);
transfer_block.copy_from_mg_add(mgdof, uu, v);
u.add(-1., uu);
deallog << u.l2_norm() << std::endl;
-
+
}
deallog.attach(logfile);
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
-
+
std::vector<double> factors;
FE_DGQ<2> q0(0);
FE_Q<2> cq1(1);
FE_RaviartThomasNodal<2> rt0(0);
FE_RaviartThomasNodal<2> rt1(1);
-
+
FESystem<2> fe0(rt1, 1, q1, 1);
FESystem<2> fe1(rt0, 2, q0, 2);
check_block(q0);
check_block(cq1);
check_block(rt1);
-
+
check_block(fe0);
check_block(fe1);
}