std::vector<std::set<unsigned int> >& boundary_indices,
const std::vector<bool>& component_mask = std::vector<bool>());
/**
+ * Maybe no longer needed.
*/
template <typename number>
extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
std::vector<std::vector<bool> > &interface_dofs,
std::vector<std::vector<bool> > &boundary_interface_dofs);
+
+ /**
+ * Does the same as the function above,
+ * but fills only the interface_dofs.
+ */
+ template <int dim, int spacedim>
+ static
+ void
+ extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
+ std::vector<std::vector<bool> > &interface_dofs);
};
/*@}*/
* matrices for each level.
*/
template <int dim, int spacedim>
- void build_matrices (const MGDoFHandler<dim,spacedim> &mg_dof);
+ void build_matrices (const MGDoFHandler<dim,spacedim> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices
+ );
virtual void prolongate (const unsigned int to_level,
VECTOR &dst,
* The data is first the global
* index, then the level index.
*/
- std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
+ std::vector<std::map<unsigned int, unsigned int> >
copy_indices;
/**
*/
std::vector<unsigned int> component_to_block_map;
- /**
- * Fill the two boolean vectors
- * #dofs_on_refinement_edge and
- * #dofs_on_refinement_boundary.
- */
- template <int dim, int spacedim>
- void find_dofs_on_refinement_edges (
- const MGDoFHandler<dim,spacedim>& mg_dof);
-
/**
* Degrees of freedom on the
* refinement edge excluding
* those on the boundary.
- *
- * @todo Clean up names
- */
- std::vector<std::vector<bool> > dofs_on_refinement_edge;
- /**
- * The degrees of freedom on the
- * the refinement edges. For each
- * level (outer vector) and each
- * dof index (inner vector), this
- * bool is true if the level
- * degree of freedom is on the
- * refinement edge towards the
- * lower level.
- *
- * @todo Clean up names
*/
- std::vector<std::vector<bool> > dofs_on_refinement_boundary;
+ std::vector<std::vector<bool> > interface_dofs;
/**
* The constraints of the global
* system.
/* --------------------- MGTransferPrebuilt -------------- */
-typedef std::vector<std::pair<unsigned int, unsigned int> >::const_iterator IT;
dst = 0;
for (unsigned int level=0;level<mg_dof_handler.get_tria().n_levels();++level)
{
+ typedef std::map<unsigned int, unsigned int>::const_iterator IT;
for (IT i= copy_indices[level].begin();
i != copy_indices[level].end();++i)
dst(i->first) = src[level](i->second);
// have fine level basis
// functions
for (unsigned int level=0;level<mg_dof_handler.get_tria().n_levels();++level)
+ {
+ typedef std::map<unsigned int, unsigned int>::const_iterator IT;
for (IT i= copy_indices[level].begin();
i != copy_indices[level].end();++i)
dst(i->first) += src[level](i->second);
+ }
}
component_to_block_map = map;
}
-
-
-template <class VECTOR>
-template <int dim, int spacedim>
-void MGTransferPrebuilt<VECTOR>::find_dofs_on_refinement_edges (
- const MGDoFHandler<dim,spacedim>& mg_dof)
-{
- for (unsigned int level = 0; level<mg_dof.get_tria().n_levels(); ++level)
- {
- std::vector<bool> tmp (mg_dof.n_dofs(level));
- dofs_on_refinement_edge.push_back (tmp);
- dofs_on_refinement_boundary.push_back (tmp);
- }
-
- const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
- const unsigned int dofs_per_face = mg_dof.get_fe().dofs_per_face;
-
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- std::vector<unsigned int> face_dof_indices (dofs_per_face);
-
- typename MGDoFHandler<dim>::cell_iterator cell = mg_dof.begin(),
- endc = mg_dof.end();
-
- for (; cell!=endc; ++cell)
- {
- std::vector<bool> cell_dofs(dofs_per_cell);
- std::vector<bool> boundary_cell_dofs(dofs_per_cell);
- const unsigned int level = cell->level();
- cell->get_mg_dof_indices (local_dof_indices);
- for (unsigned int face_nr=0;
- face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
- {
- typename DoFHandler<dim>::face_iterator face = cell->face(face_nr);
- if(!cell->at_boundary(face_nr))
- {
- //interior face
- typename MGDoFHandler<dim>::cell_iterator neighbor
- = cell->neighbor(face_nr);
- // Do refinement face
- // from the coarse side
- if (neighbor->level() < cell->level())
- {
- for(unsigned int j=0; j<dofs_per_face; ++j)
- cell_dofs[mg_dof.get_fe().face_to_cell_index(j,face_nr)] = true;
- }
- }
- //boundary face
- else
- for(unsigned int j=0; j<dofs_per_face; ++j)
- boundary_cell_dofs[mg_dof.get_fe().face_to_cell_index(j,face_nr)] = true;
- }//faces
- for(unsigned int i=0; i<dofs_per_cell; ++i)
- {
- if(cell_dofs[i])
- {
- dofs_on_refinement_edge[level][local_dof_indices[i]] = true;
- dofs_on_refinement_boundary[level][local_dof_indices[i]] = true;
- }
- if(boundary_cell_dofs[i])
- dofs_on_refinement_boundary[level][local_dof_indices[i]] = true;
- }
- }
-}
-
template <class VECTOR>
unsigned int
MGTransferPrebuilt<VECTOR>::memory_consumption () const
#include <base/config.h>
#include <lac/block_vector.h>
+#include <lac/constraint_matrix.h>
#ifdef DEAL_PREFER_MATRIX_EZ
# include <lac/sparse_matrix_ez.h>
# include <lac/block_sparse_matrix_ez.h>
* Start index of each component.
*/
std::vector<unsigned int> component_start;
-
+
/**
* Start index of each component on
* all levels.
*/
std::vector<std::vector<unsigned int> > mg_component_start;
-
+
/**
* Call build_matrices()
* function first.
std::vector<std_cxx1x::shared_ptr<BlockSparseMatrix<double> > > prolongation_matrices;
/**
- * Unused now, but intended to
- * hold the mapping for the
+ * Holds the mapping for the
* <tt>copy_to/from_mg</tt>-functions.
+ * The data is first the global
+ * index, then the level index.
*/
std::vector<std::map<unsigned int, unsigned int> >
copy_to_and_from_indices;
+
+ /**
+ * Store the boundary_indices.
+ * These are needed for the
+ * boundary values in the
+ * restriction matrix.
+ */
+ std::vector<std::set<unsigned int> > boundary_indices;
};
//TODO:[GK] Update documentation for copy_* functions
private MGTransferComponentBase
{
public:
+ /**
+ * Constructor without constraint
+ * matrices. Use this constructor
+ * only with discontinuous finite
+ * elements or with no local
+ * refinement.
+ */
+ MGTransferSelect ();
+
+ /**
+ * Constructor with constraint matrices.
+ */
+ MGTransferSelect (const ConstraintMatrix& constraints);
+
/**
* Destructor.
*/
* DoFRenumbering::component_wise). It
* also affects the behavior of
* the <tt>selected</tt> argument
+ *
+ * @arg boundary_indices: holds the
+ * boundary indices on each level.
*/
template <int dim, int spacedim>
void build_matrices (const DoFHandler<dim,spacedim> &dof,
const std::vector<unsigned int>& target_component
= std::vector<unsigned int>(),
const std::vector<unsigned int>& mg_target_component
- = std::vector<unsigned int>());
+ = std::vector<unsigned int>(),
+ const std::vector<std::set<unsigned int> >& boundary_indices
+ = std::vector<std::set<unsigned int> >()
+ );
/**
* Change selected
void
do_copy_from_mg (const MGDoFHandler<dim,spacedim> &mg_dof,
OutVector &dst,
- const MGLevelObject<Vector<number> > &src,
- const unsigned int offset) const;
+ const MGLevelObject<Vector<number> > &src) const;
/**
* Implementation of the public
void
do_copy_from_mg_add (const MGDoFHandler<dim,spacedim> &mg_dof,
OutVector &dst,
- const MGLevelObject<Vector<number> > &src,
- const unsigned int offset) const;
+ const MGLevelObject<Vector<number> > &src) const;
/**
* Actual implementation of
void
do_copy_to_mg (const MGDoFHandler<dim,spacedim> &mg_dof,
MGLevelObject<Vector<number> > &dst,
- const InVector &src,
- const unsigned int offset) const;
+ const InVector &src) const;
/**
* Selected component of global vector.
*/
* Selected component inside multigrid.
*/
unsigned int mg_selected_component;
+
+ /**
+ * The degrees of freedom on the
+ * the refinement edges. For each
+ * level (outer vector) and each
+ * dof index (inner vector), this
+ * bool is true if the level
+ * degree of freedom is on the
+ * refinement edge towards the
+ * lower level excluding boundary dofs.
+ */
+ std::vector<std::vector<bool> > interface_dofs;
+
+ /**
+ * The constraints of the global
+ * system.
+ */
+ public:
+ SmartPointer<const ConstraintMatrix> constraints;
};
/*@}*/
#include <multigrid/mg_dof_accessor.h>
#include <multigrid/mg_tools.h>
#include <multigrid/mg_transfer_component.h>
+#include <numerics/data_out.h>
#include <algorithm>
+#include <sstream>
+#include <fstream>
DEAL_II_NAMESPACE_OPEN
MGLevelObject<Vector<number> > &dst,
const BlockVector<number2> &src) const
{
- do_copy_to_mg (mg_dof_handler, dst, src, 0);
+ do_copy_to_mg (mg_dof_handler, dst, src);
}
MGLevelObject<Vector<number> > &dst,
const Vector<number2> &src) const
{
- do_copy_to_mg (mg_dof_handler, dst, src, component_start[selected_component]);
+ do_copy_to_mg (mg_dof_handler, dst, src);
}
BlockVector<number2>& dst,
const MGLevelObject<Vector<number> >& src) const
{
- do_copy_from_mg (mg_dof_handler, dst, src, 0);
+ do_copy_from_mg (mg_dof_handler, dst, src);
+ if (constraints != 0)
+ constraints->condense(dst);
}
Vector<number2>& dst,
const MGLevelObject<Vector<number> >& src) const
{
- do_copy_from_mg (mg_dof_handler, dst, src,
- component_start[selected_component]);
+ do_copy_from_mg (mg_dof_handler, dst, src);
+ if (constraints != 0)
+ {
+ //If we were given constraints
+ //apply them to the dst that goes
+ //back now to the linear solver.
+ //Since constraints are globally
+ //defined create a global vector here
+ //and copy dst to the right component,
+ //apply the constraints then and copy
+ //the block back to dst.
+ const unsigned int n_blocks =
+ *std::max_element(target_component.begin(), target_component.end()) + 1;
+ std::vector<unsigned int> dofs_per_block (n_blocks);
+ DoFTools::count_dofs_per_block (mg_dof_handler, dofs_per_block, target_component);
+ BlockVector<number> tmp;
+ tmp.reinit(n_blocks);
+ for(unsigned int b=0; b<n_blocks; ++b)
+ tmp.block(b).reinit(dofs_per_block[b]);
+ tmp.collect_sizes ();
+ tmp.block(target_component[selected_component]) = dst;
+ constraints->condense(tmp);
+ dst = tmp.block(target_component[selected_component]);
+ }
}
BlockVector<number2>& dst,
const MGLevelObject<Vector<number> >& src) const
{
- do_copy_from_mg_add (mg_dof_handler, dst, src, 0);
+ do_copy_from_mg_add (mg_dof_handler, dst, src);
}
Vector<number2>& dst,
const MGLevelObject<Vector<number> >& src) const
{
- do_copy_from_mg_add (mg_dof_handler, dst, src,
- component_start[selected_component]);
+ do_copy_from_mg_add (mg_dof_handler, dst, src);
}
MGTransferSelect<number>::do_copy_from_mg (
const MGDoFHandler<dim,spacedim> &mg_dof_handler,
OutVector &dst,
- const MGLevelObject<Vector<number> > &src,
- const unsigned int offset) const
+ const MGLevelObject<Vector<number> > &src) const
{
- const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+// const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- std::vector<unsigned int> global_dof_indices (dofs_per_cell);
- std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+// const unsigned int dofs_per_cell = fe.dofs_per_cell;
+// std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+// std::vector<unsigned int> level_dof_indices (dofs_per_cell);
typename MGDoFHandler<dim,spacedim>::active_cell_iterator
level_cell = mg_dof_handler.begin_active();
// Note that the level is
// monotonuosly increasing
+ dst = 0;
for (; level_cell != endc; ++level_cell)
- {
- const unsigned int level = level_cell->level();
-
- // get the dof numbers of
- // this cell for the global
- // and the level-wise
- // numbering
- level_cell->get_dof_indices (global_dof_indices);
- level_cell->get_mg_dof_indices(level_dof_indices);
-
- // copy level-wise data to
- // global vector
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const unsigned int component
- = mg_target_component[fe.system_to_component_index(i).first];
- if (mg_selected[component])
- {
- const unsigned int level_start
- = mg_component_start[level][component];
- dst(global_dof_indices[i] - offset)
- = src[level](level_dof_indices[i]-level_start);
- }
- }
- }
+ {
+ const unsigned int level = level_cell->level();
+ typedef std::map<unsigned int, unsigned int>::const_iterator IT;
+ for (IT i=copy_to_and_from_indices[level].begin();
+ i != copy_to_and_from_indices[level].end(); ++i)
+ dst(i->first) = src[level](i->second);
+ }
}
-
template <typename number>
template <int dim, class OutVector, int spacedim>
void
MGTransferSelect<number>::do_copy_from_mg_add (
const MGDoFHandler<dim,spacedim> &mg_dof_handler,
OutVector &dst,
- const MGLevelObject<Vector<number> > &src,
- const unsigned int offset) const
+ const MGLevelObject<Vector<number> > &src) const
{
const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
// Note that the level is
// monotonuosly increasing
+ dst = 0;
for (; level_cell != endc; ++level_cell)
{
const unsigned int level = level_cell->level();
-
- // get the dof numbers of
- // this cell for the global
- // and the level-wise
- // numbering
- level_cell->get_dof_indices (global_dof_indices);
- level_cell->get_mg_dof_indices(level_dof_indices);
- // copy level-wise data to
- // global vector
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const unsigned int component
- = mg_target_component[fe.system_to_component_index(i).first];
- if (mg_selected[component])
- {
- const unsigned int level_start
- = mg_component_start[level][component];
- dst(global_dof_indices[i] - offset)
- += src[level](level_dof_indices[i] - level_start);
- }
- }
+ typedef std::map<unsigned int, unsigned int>::const_iterator IT;
+ for (IT i=copy_to_and_from_indices[level].begin();
+ i != copy_to_and_from_indices[level].end(); ++i)
+ dst(i->first) += src[level](i->second);
}
}
#include <base/logstream.h>
-#include "iostream"
+#include <iostream>
DEAL_II_NAMESPACE_OPEN
deallog << "Smoothing on level " << level << std::endl;
// smoothing of the residual by
// modifying s
+// defect[level].print(std::cout, 2,false);
+// std::cout<<std::endl;
pre_smooth->smooth(level, solution[level], defect[level]);
-
+// solution[level].print(std::cout, 2,false);
+
if (debug>2)
deallog << "Solution norm " << solution[level].l2_norm()
<< std::endl;
deallog << "Residual on level " << level << std::endl;
// t = A*solution[level]
matrix->vmult(level, t[level], solution[level]);
+// std::cout<<std::endl;
+// t[level].print(std::cout, 2,false);
// make t rhs of lower level The
// non-refined parts of the
{
t[l-1] = 0.;
if (l==level && edge_out != 0)
- edge_out->vmult_add(level, t[level], solution[level]);
-
+ {
+ edge_out->vmult_add(level, t[level], solution[level]);
+ deallog << "Norm t[" << level << "] " << t[level].l2_norm() << std::endl;
+ }
+
if (l==level && edge_down != 0)
edge_down->vmult(level, t[level-1], solution[level]);
transfer->restrict_and_add (l, t[l-1], t[l]);
+ deallog << "restrict t[" << l-1 << "] " << t[l-1].l2_norm() << std::endl;
defect[l-1] -= t[l-1];
+ deallog << "defect d[" << l-1 << "] " << defect[l-1].l2_norm() << std::endl;
}
// do recursion
}
if (debug>2)
- deallog << "V-cycle Defect norm " << defect[level].l2_norm()
+ deallog << "V-cycle Defect norm hier " << defect[level].l2_norm()
<< std::endl;
if (debug>1)
deallog << "Smoothing on level " << level << std::endl;
// post-smoothing
+
+// std::cout<<std::endl;
+// defect[level].print(std::cout, 2,false);
post_smooth->smooth(level, solution[level], defect[level]);
-
+// solution[level].print(std::cout, 2,false);
+// std::cout<<std::endl;
+
if (debug>2)
deallog << "Solution norm " << solution[level].l2_norm()
<< std::endl;