#include <dofs/dof_tools.h>
#include <fe/fe.h>
-#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
std::vector<std::vector<bool> >
dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
std::vector<std::vector<bool> >
- block_select (n_blocks, std::vector<bool>(n_blocks, false));
+ block_select;
+ block_select.resize(n_blocks);
+ for (int iii=0;iii<n_blocks;++iii)
+ block_select[iii].resize(n_blocks, false);
Threads::TaskGroup<> tasks;
for (unsigned int i=0; i<n_blocks; ++i)
{
+template <int dim, int spacedim>
+void
+MGTools::
+extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
+ std::vector<std::vector<bool> > &interface_dofs)
+{
+ Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_levels(),
+ ExcDimensionMismatch (interface_dofs.size(),
+ mg_dof_handler.get_tria().n_levels()));
+
+ for (unsigned int l=0; l<mg_dof_handler.get_tria().n_levels(); ++l)
+ {
+ Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
+ ExcDimensionMismatch (interface_dofs[l].size(),
+ mg_dof_handler.n_dofs(l)));
+
+ std::fill (interface_dofs[l].begin(),
+ interface_dofs[l].end(),
+ false);
+ }
+
+ const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<bool> cell_dofs(dofs_per_cell, false);
+
+ typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+ endc = mg_dof_handler.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ std::fill (cell_dofs.begin(), cell_dofs.end(), false);
+
+ for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+ {
+ const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
+ if (!face->at_boundary())
+ {
+ //interior face
+ const typename MGDoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(face_nr);
+
+ if (neighbor->level() < cell->level())
+ {
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+
+ }
+ }
+ }
+
+ const unsigned int level = cell->level();
+ cell->get_mg_dof_indices (local_dof_indices);
+
+ for(unsigned int i=0; i<dofs_per_cell; ++i)
+ if (cell_dofs[i])
+ interface_dofs[level][local_dof_indices[i]] = true;
+ }
+}
template <int dim, int spacedim>
for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
if(cell->at_boundary(face_nr))
for(unsigned int j=0; j<dofs_per_face; ++j)
- if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true)
+ if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) //is this necessary?
boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
std::vector<std::vector<bool> > &interface_dofs,
std::vector<std::vector<bool> > &boundary_interface_dofs);
+template
+void
+MGTools::
+extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
+ std::vector<std::vector<bool> > &interface_dofs);
#endif
DEAL_II_NAMESPACE_CLOSE
//---------------------------------------------------------------------------
#include <base/logstream.h>
+#include <base/function.h>
#include <lac/vector.h>
#include <lac/block_vector.h>
+#include <lac/block_indices.h>
#include <lac/sparse_matrix.h>
#include <lac/block_sparse_matrix.h>
#include <grid/tria.h>
#include <algorithm>
#include <numeric>
-
+#include <iostream>
DEAL_II_NAMESPACE_OPEN
new_dofs(mg_dof.get_tria().n_levels(),
std::vector<unsigned int>(target_component.size()));
std::swap(ndofs, new_dofs);
- MGTools::count_dofs_per_component (mg_dof, ndofs,
- true, target_component);
+ MGTools::count_dofs_per_block (mg_dof, ndofs, target_component);
}
for (unsigned int level=v.get_minlevel();
// Compute the number of blocks needed
#ifdef DEBUG
- const unsigned int n_selected
- = std::accumulate(selected.begin(),
- selected.end(),
- 0U);
- Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1));
+// const unsigned int n_selected
+// = std::accumulate(selected.begin(),
+// selected.end(),
+// 0U);
+// Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1));
#endif
unsigned int selected_block = 0;
- while (!selected[selected_block])
- ++selected_block;
-
+ for (unsigned int i=0; i<target_component.size(); ++i)
+ if(selected[i])
+ selected_block = target_component[i];
+
if (ndofs.size() == 0)
{
std::vector<std::vector<unsigned int> >
new_dofs(mg_dof.get_tria().n_levels(),
std::vector<unsigned int>(target_component.size()));
std::swap(ndofs, new_dofs);
- MGTools::count_dofs_per_component (mg_dof, ndofs,
- true, target_component);
+ MGTools::count_dofs_per_block (mg_dof, ndofs,
+ target_component);
}
for (unsigned int level=v.get_minlevel();
MGTransferSelect<number>::do_copy_to_mg (
const MGDoFHandler<dim,spacedim>& mg_dof_handler,
MGLevelObject<Vector<number> >& dst,
- const InVector& src,
- const unsigned int offset) const
+ const InVector& src) const
{
- const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
- // set the elements of the vectors
- // on all levels to zero
- unsigned int minlevel = dst.get_minlevel();
- unsigned int maxlevel = dst.get_maxlevel();
-
dst=0;
Assert(sizes.size()==mg_dof_handler.get_tria().n_levels(),
reinit_vector_by_components(mg_dof_handler, dst, mg_selected,
mg_target_component, sizes);
- std::vector<unsigned int> global_dof_indices (dofs_per_cell);
- std::vector<unsigned int> level_dof_indices (dofs_per_cell);
-
- // Build a vector of the selected
- // indices, since traversing all
- // indices on each cell is too
- // slow.
- std::vector<unsigned int> selected_indices;
- selected_indices.reserve(dofs_per_cell);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (mg_selected[mg_target_component[fe.system_to_component_index(i).first]])
- selected_indices.push_back(i);
- unsigned int selected_component
- = mg_target_component[fe.system_to_component_index(selected_indices[0]).first];
-
// traverse the grid top-down
// (i.e. starting with the most
// refined grid). this way, we can
// the respective vector on the
// next finer level, which we then
// already have built.
- for (int level=maxlevel; level>=static_cast<signed int>(minlevel); --level)
- {
- typename MGDoFHandler<dim,spacedim>::active_cell_iterator
- level_cell = mg_dof_handler.begin_active(level);
- const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
- level_end = mg_dof_handler.end_active(level);
- // Compute coarse level right hand side
- // by restricting from fine level.
- for (; level_cell!=level_end; ++level_cell)
- {
- // 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);
-
- // transfer the global
- // defect in the vector
- // into the level-wise one
- const unsigned int level_start
- = mg_component_start[level][selected_component];
- const typename std::vector<unsigned int>::const_iterator
- end = selected_indices.end();
-
- for (typename std::vector<unsigned int>::const_iterator
- i=selected_indices.begin();
- i != end; ++i)
- {
- dst[level](level_dof_indices[*i] - level_start)
- = src(global_dof_indices[*i] - offset);
- }
- }
+ bool first = true;
+ for (unsigned int level=mg_dof_handler.get_tria().n_levels(); level!=0;)
+ {
+ --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[level](i->second) = src(i->first);
// for that part of the level
// which is further refined:
// get the defect by
// restriction of the defect on
// one level higher
- if (static_cast<unsigned int>(level) < maxlevel)
- {
- restrict_and_add (level+1, dst[level], dst[level+1]);
- }
+ if(!first)
+ restrict_and_add (level+1, dst[level], dst[level+1]);
+ first = false;
}
}
// Compute the lengths of all blocks
sizes.resize(n_levels);
- MGTools::count_dofs_per_component(mg_dof, sizes, true, mg_target_component);
+ for(unsigned int l=0; l<n_levels; ++l)
+ sizes[l].resize(n_components);
+
+ MGTools::count_dofs_per_block(mg_dof, sizes, mg_target_component);
// Fill some index vectors
// for later use.
}
component_start.resize(*std::max_element (target_component.begin(),
- target_component.end()) + 1);
+ target_component.end()) + 1);
DoFTools::
- count_dofs_per_component (static_cast<const DoFHandler<dim,spacedim>&>(mg_dof),
- component_start, true, target_component);
+ count_dofs_per_block (mg_dof, component_start, target_component);
unsigned int k=0;
for (unsigned int i=0;i<component_start.size();++i)
= fe.system_to_component_index(i).first;
const unsigned int jcomp
= fe.system_to_component_index(j).first;
- if ((icomp==jcomp) && mg_selected[mg_target_component[icomp]])
+ if ((icomp==jcomp) && mg_selected[icomp])
prolongation_sparsities[level]->add(dof_indices_child[i],
dof_indices_parent[j]);
};
{
const unsigned int icomp = fe.system_to_component_index(i).first;
const unsigned int jcomp = fe.system_to_component_index(j).first;
- if ((icomp==jcomp) && mg_selected[mg_target_component[icomp]])
+ if ((icomp==jcomp) && mg_selected[icomp])
prolongation_matrices[level]->set(dof_indices_child[i],
dof_indices_parent[j],
prolongation(i,j));
}
}
}
+ // impose boundary conditions
+ // but only in the column of
+ // the prolongation matrix
+ //TODO: this way is not very efficient
+
+ if(boundary_indices.size() != 0)
+ {
+ std::vector<std::vector<unsigned int> >
+ dofs_per_component(mg_dof.get_tria().n_levels(),
+ std::vector<unsigned int>(n_components));
+
+ MGTools::count_dofs_per_block (mg_dof, dofs_per_component, mg_target_component);
+ for (unsigned int level=0; level<n_levels-1; ++level)
+ {
+ if (boundary_indices[level].size() == 0)
+ continue;
+
+ for (unsigned int iblock=0; iblock<n_components; ++iblock)
+ for (unsigned int jblock=0; jblock<n_components; ++jblock)
+ if (iblock==jblock)
+ {
+ const unsigned int n_dofs = prolongation_matrices[level]->block(iblock,jblock).m();
+ for (unsigned int i=0; i<n_dofs; ++i)
+ {
+ SparseMatrix<double>::iterator anfang = prolongation_matrices[level]->block(iblock,jblock).begin(i),
+ ende = prolongation_matrices[level]->block(iblock,jblock).end(i);
+ for(; anfang != ende; ++anfang)
+ {
+ const unsigned int column_number = anfang->column();
+
+ //convert global indices into local ones
+ const BlockIndices block_indices_coarse (dofs_per_component[level]);
+ const unsigned int global_j = block_indices_coarse.local_to_global(iblock, column_number);
+
+ std::set<unsigned int>::const_iterator found_dof =
+ boundary_indices[level].find(global_j);
+
+ const bool is_boundary_index =
+ (found_dof != boundary_indices[level].end());
+
+ if(is_boundary_index)
+ {
+ prolongation_matrices[level]->block(iblock,jblock)
+ .set(i,column_number,0);
+ }
+ }
+ }
+ }
+ }
+ }
}
unsigned int select,
unsigned int mg_select,
const std::vector<unsigned int>& t_component,
- const std::vector<unsigned int>& mg_t_component)
+ const std::vector<unsigned int>& mg_t_component,
+ const std::vector<std::set<unsigned int> >& bdry_indices)
{
const FiniteElement<dim>& fe = mg_dof.get_fe();
unsigned int ncomp = mg_dof.get_fe().n_components();
target_component = t_component;
mg_target_component = mg_t_component;
+ boundary_indices = bdry_indices;
selected_component = select;
mg_selected_component = mg_select;
selected.resize(ncomp, false);
- selected[select] = true;
+ for(unsigned int c=0; c<ncomp; ++c)
+ if(t_component[c] == selected_component)
+ selected[c] = true;
+
mg_selected.resize(ncomp, false);
- mg_selected[mg_select] = true;
+ for(unsigned int c=0; c<ncomp; ++c)
+ if(mg_t_component[c] == mg_selected_component)
+ mg_selected[c] = true;
// If components are renumbered,
// find the first original
// component corresponding to the
break;
}
}
+
MGTransferComponentBase::build_matrices (dof, mg_dof);
+ interface_dofs.resize(mg_dof.get_tria().n_levels());
+ for(unsigned int l=0; l<mg_dof.get_tria().n_levels(); ++l)
+ interface_dofs[l].resize(mg_dof.n_dofs(l));
+ MGTools::extract_inner_interface_dofs(mg_dof, interface_dofs);
+
std::vector<unsigned int> global_dof_indices (fe.dofs_per_cell);
std::vector<unsigned int> level_dof_indices (fe.dofs_per_cell);
for (int level=dof.get_tria().n_levels()-1; level>=0; --level)
{
+ copy_to_and_from_indices[level].clear();
typename MGDoFHandler<dim,spacedim>::active_cell_iterator
level_cell = mg_dof.begin_active(level);
const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
{
const unsigned int component
- = mg_target_component[fe.system_to_component_index(i).first];
- if (mg_selected[component])
+ = fe.system_to_component_index(i).first;
+ if (mg_selected[component] &&
+ !interface_dofs[level][level_dof_indices[i]])
{
const unsigned int level_start
- = mg_component_start[level][component];
+ = mg_component_start[level][mg_target_component[component]];
copy_to_and_from_indices[level].insert(
std::make_pair(global_dof_indices[i]
- - component_start[selected_component],
+ - component_start[target_component[selected_component]],
level_dof_indices[i]-level_start));
}
}
const MGDoFHandler<deal_II_dimension> &,
unsigned int, unsigned int,
const std::vector<unsigned int>&,
- const std::vector<unsigned int>&);
+ const std::vector<unsigned int>&,
+ const std::vector<std::set<unsigned int> >&);
template
void MGTransferSelect<double>::build_matrices<deal_II_dimension>
const MGDoFHandler<deal_II_dimension> &,
unsigned int, unsigned int,
const std::vector<unsigned int>&,
- const std::vector<unsigned int>&);
+ const std::vector<unsigned int>&,
+ const std::vector<std::set<unsigned int> >&);
template void
MGTransferSelect<float>::copy_to_mg (
--level;
VECTOR& dst_level = dst[level];
- typedef std::vector<std::pair<unsigned int, unsigned int> >::const_iterator IT;
+ typedef std::map<unsigned int, unsigned int>::const_iterator IT;
for (IT i= copy_indices[level].begin();
i != copy_indices[level].end();++i)
dst_level(i->second) = src(i->first);
template <typename VECTOR>
template <int dim, int spacedim>
void MGTransferPrebuilt<VECTOR>::build_matrices (
- const MGDoFHandler<dim,spacedim> &mg_dof)
+ const MGDoFHandler<dim,spacedim> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices
+ )
{
//start building the matrices here
const unsigned int n_levels = mg_dof.get_tria().n_levels();
// impose boundary conditions
// but only in the column of
// the prolongation matrix
- std::vector<std::set<unsigned int> > boundary_indices (n_levels);
- typename FunctionMap<dim>::type boundary;
- ZeroFunction<dim> homogeneous_dirichlet_bc (1);
- boundary[0] = &homogeneous_dirichlet_bc;
- MGTools::make_boundary_list(mg_dof, boundary, boundary_indices);
-
- std::vector<unsigned int> constrain_indices;
- for (int level=n_levels-2; level>=0; --level)
+
+ if(boundary_indices.size() != 0)
+ {
+ std::vector<unsigned int> constrain_indices;
+ for (int level=n_levels-2; level>=0; --level)
+// for (unsigned int level=0; level<n_levels-1; ++level)
{
if (boundary_indices[level].size() == 0)
- continue;
+ continue;
- // need to delete all the columns in the
- // matrix that are on the boundary. to fix
- // this, create an array as long as there are
- // matrix columns, and find which columns we
- // need to filter away.
+ // need to delete all the columns in the
+ // matrix that are on the boundary. to fix
+ // this, create an array as long as there are
+ // matrix columns, and find which columns we
+ // need to filter away.
constrain_indices.resize (0);
constrain_indices.resize (prolongation_matrices[level]->n(), 0);
std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
- endd = boundary_indices[level].end();
+ endd = boundary_indices[level].end();
for (; dof != endd; ++dof)
- constrain_indices[*dof] = 1;
+ constrain_indices[*dof] = 1;
const unsigned int n_dofs = prolongation_matrices[level]->m();
for (unsigned int i=0; i<n_dofs; ++i)
- {
- SparseMatrix<double>::iterator start_row = prolongation_matrices[level]->begin(i),
- end_row = prolongation_matrices[level]->end(i);
- for(; start_row != end_row; ++start_row)
- {
- if (constrain_indices[start_row->column()] == 1)
- start_row->value() = 0;
- }
- }
+ {
+ SparseMatrix<double>::iterator start_row = prolongation_matrices[level]->begin(i),
+ end_row = prolongation_matrices[level]->end(i);
+ for(; start_row != end_row; ++start_row)
+ {
+ std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
+ endd = boundary_indices[level].end();
+ for (; dof != endd; ++dof)
+ {
+ const unsigned int column_number = start_row->column();
+ const unsigned int dof_number = *dof;
+ if(dof_number == column_number)
+ prolongation_matrices[level]->set(i,dof_number,0);
+ }
+ }
+ }
}
+ }
// to find the indices that describe the
// relation between global dofs and local
// the std::vector<std::pair<uint,uint> > that
// only contains the active dofs on the
// levels.
-
- find_dofs_on_refinement_edges (mg_dof);
+ interface_dofs.resize(mg_dof.get_tria().n_levels());
+ for(unsigned int l=0; l<mg_dof.get_tria().n_levels(); ++l)
+ interface_dofs[l].resize(mg_dof.n_dofs(l));
+ MGTools::extract_inner_interface_dofs (mg_dof, interface_dofs);
copy_indices.resize(n_levels);
- std::vector<unsigned int> temp_copy_indices;
+// std::vector<unsigned int> temp_copy_indices;
std::vector<unsigned int> global_dof_indices (dofs_per_cell);
std::vector<unsigned int> level_dof_indices (dofs_per_cell);
for (int level=mg_dof.get_tria().n_levels()-1; level>=0; --level)
const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
level_end = mg_dof.end_active(level);
- temp_copy_indices.resize (0);
- temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int);
+// temp_copy_indices.resize (0);
+// temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int);
// Compute coarse level right hand side
// by restricting from fine level.
level_cell->get_mg_dof_indices (level_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
- if(!dofs_on_refinement_edge[level][level_dof_indices[i]])
- temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+ {
+ if(!interface_dofs[level][level_dof_indices[i]])
+ copy_indices[level].insert(
+ std::make_pair(global_dof_indices[i], level_dof_indices[i]));
+// temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+// if(!interface_dofs[level][level_dof_indices[i]])
+// temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+ }
}
// now all the active dofs got a valid entry,
// copy_indices object. Then, insert the pairs
// of global index and level index into
// copy_indices.
- const unsigned int n_active_dofs =
- std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(),
- std::bind2nd(std::not_equal_to<unsigned int>(),
- numbers::invalid_unsigned_int));
- copy_indices[level].resize (n_active_dofs);
- unsigned int counter = 0;
- for (unsigned int i=0; i<temp_copy_indices.size(); ++i)
- if (temp_copy_indices[i] != numbers::invalid_unsigned_int)
- copy_indices[level][counter++] =
- std::make_pair<unsigned int,unsigned int> (temp_copy_indices[i], i);
- Assert (counter == n_active_dofs, ExcInternalError());
+// const unsigned int n_active_dofs =
+// std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(),
+// std::bind2nd(std::not_equal_to<unsigned int>(),
+// numbers::invalid_unsigned_int));
+// copy_indices[level].resize (n_active_dofs);
+// unsigned int counter = 0;
+// for (unsigned int i=0; i<temp_copy_indices.size(); ++i)
+// if (temp_copy_indices[i] != numbers::invalid_unsigned_int)
+// copy_indices[level][counter++] =
+// std::make_pair<unsigned int,unsigned int> (temp_copy_indices[i], i);
+// Assert (counter == n_active_dofs, ExcInternalError());
}
}
template
void MGTransferPrebuilt<Vector<float> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
template
void MGTransferPrebuilt<Vector<double> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
template
void MGTransferPrebuilt<BlockVector<float> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
template
void MGTransferPrebuilt<BlockVector<double> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
template void
MGTransferPrebuilt<Vector<float> >::copy_to_mg (