#include <iostream>
#include <fstream>
#include <sstream>
-#include <boost/archive/text_oarchive.hpp>
-#include <boost/archive/text_iarchive.hpp>
// The last step is as in all
// previous programs:
using namespace dealii;
-
-
-
-
-
// @sect3{The <code>LaplaceProblem</code> class template}
// This main class is basically the same
const unsigned int) const
{
if (p.square() < 0.5*0.5)
- return 20;
+ return 5;
else
return 1;
}
deallog.depth_console(0);
}
- std::string id_to_string(const CellId &id)
- {
- std::ostringstream ss;
- ss << id;
- return ss.str();
- }
-
// @sect4{LaplaceProblem::setup_system}
// freedom at once can be done using
// ConstraintMatrix::add_lines():
std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_global_levels());
- std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_global_levels());
+ ConstraintMatrix empty_constraints;
for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
{
IndexSet dofset;
boundary_constraints[level].close ();
-
- boundary_interface_constraints[level].reinit(dofset);
- // compute indices that are refinement edge and boundary (aka refinement_edge_boundary_indices):
- IndexSet idxset = mg_constrained_dofs.get_refinement_edge_indices(level) & mg_constrained_dofs.get_boundary_indices(level);
-
- boundary_interface_constraints[level]
- .add_lines (idxset);
- boundary_interface_constraints[level].close ();
}
// Now that we're done with most of our preliminaries, let's start
const IndexSet &interface_dofs_on_level
= mg_constrained_dofs.get_refinement_edge_indices(cell->level());
-
+ const unsigned int lvl = cell->level();
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- /** old HEAD:
- if (!mg_constrained_dofs.at_refinement_edge(cell->level(),local_dof_indices[i])
- || mg_constrained_dofs.at_refinement_edge(cell->level(),local_dof_indices[j])) */
- if ( !(interface_dofs_on_level.is_element(local_dof_indices[i])==true &&
- interface_dofs_on_level.is_element(local_dof_indices[j])==false))
- cell_matrix(i,j) = 0;
-
- boundary_interface_constraints[cell->level()]
+ if (interface_dofs_on_level.is_element(local_dof_indices[i]) // at_refinement_edge(i)
+ &&
+ !interface_dofs_on_level.is_element(local_dof_indices[j]) // !at_refinement_edge(j)
+ &&
+ (
+ (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+ ) // ( !boundary(i) && !boundary(j) )
+ ||
+ (
+ mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ local_dof_indices[i]==local_dof_indices[j]
+ ) // ( boundary(i) && boundary(j) && i==j )
+ )
+ )
+ {
+ // do nothing, so add entries to interface matrix
+ }
+ else
+ {
+ cell_matrix(i,j) = 0;
+ }
+
+
+ empty_constraints
.distribute_local_to_global (cell_matrix,
local_dof_indices,
mg_interface_matrices[cell->level()]);
SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false);
SolverGMRES<vector_t> cg (solver_control);
- //solution = 0;
-
if (false)
{
// code to optionally compare to Trilinos ML
template <int dim>
void LaplaceProblem<dim>::refine_grid ()
{
-
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
TrilinosWrappers::MPI::Vector temp_solution;
estimated_error_per_cell,
0.3, 0.0);
-
- triangulation.prepare_coarsening_and_refinement ();
triangulation.execute_coarsening_and_refinement ();
}
? ")" : ", ");
deallog << std::endl;
-
assemble_system ();
assemble_multigrid ();
= mg_constrained_dofs.get_refinement_edge_boundary_indices ();
std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
- std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
+ ConstraintMatrix empty_constraints;
for (unsigned int level=0; level<triangulation.n_levels(); ++level)
{
boundary_constraints[level].add_lines (interface_dofs[level]);
boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
boundary_constraints[level].close ();
-
- boundary_interface_constraints[level]
- .add_lines (boundary_interface_dofs[level]);
- boundary_interface_constraints[level].close ();
}
typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
.distribute_local_to_global (cell_matrix,
local_dof_indices,
mg_matrices[cell->level()]);
+ const unsigned int lvl = cell->level();
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ( !(interface_dofs[cell->level()][local_dof_indices[i]]==true &&
- interface_dofs[cell->level()][local_dof_indices[j]]==false))
- cell_matrix(i,j) = 0;
-
- boundary_interface_constraints[cell->level()]
+ if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i)
+ &&
+ !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j)
+ &&
+ (
+ (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+ ) // ( !boundary(i) && !boundary(j) )
+ ||
+ (
+ mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ local_dof_indices[i]==local_dof_indices[j]
+ ) // ( boundary(i) && boundary(j) && i==j )
+ )
+ )
+ {
+ // do nothing, so add entries to interface matrix
+ }
+ else
+ {
+ cell_matrix(i,j) = 0;
+ }
+
+ empty_constraints
.distribute_local_to_global (cell_matrix,
local_dof_indices,
mg_interface_matrices[cell->level()]);
private:
void setup_system ();
void assemble_system ();
- void assemble_multigrid (const bool &use_mw);
- void solve ();
+ void assemble_multigrid (bool use_mw);
+ void solve (bool use_mw);
void refine_grid (const std::string &reftype);
void output_results (const unsigned int cycle) const;
template <int dim>
-void LaplaceProblem<dim>::assemble_multigrid (const bool &use_mw)
+void LaplaceProblem<dim>::assemble_multigrid (bool use_mw)
{
+ deallog << "assemble_multigrid " << (use_mw ? "(mesh_worker)":"") << std::endl;
+
if (use_mw == true)
{
mg_matrices = 0.;
+ mg_interface_in = 0.;
+ mg_interface_out = 0.;
MappingQ1<dim> mapping;
MeshWorker::IntegrationInfoBox<dim> info_box;
std::vector<std::vector<bool> > interface_dofs
= mg_constrained_dofs.get_refinement_edge_indices ();
- std::vector<std::vector<bool> > boundary_interface_dofs
- = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
- std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
+ ConstraintMatrix empty_constraints;
+
for (unsigned int level=0; level<triangulation.n_levels(); ++level)
{
boundary_constraints[level].add_lines (interface_dofs[level]);
boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
boundary_constraints[level].close ();
-
- boundary_interface_constraints[level]
- .add_lines (boundary_interface_dofs[level]);
- boundary_interface_constraints[level].close ();
}
typename DoFHandler<dim>::level_cell_iterator cell = mg_dof_handler.begin_mg(),
local_dof_indices,
mg_matrices[cell->level()]);
+ const unsigned int lvl = cell->level();
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ( !(interface_dofs[cell->level()][local_dof_indices[i]]==true &&
- interface_dofs[cell->level()][local_dof_indices[j]]==false))
- cell_matrix(i,j) = 0;
+ if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i)
+ &&
+ !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j)
+ &&
+ (
+ (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+ ) // ( !boundary(i) && !boundary(j) )
+ ||
+ (
+ mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ local_dof_indices[i]==local_dof_indices[j]
+ ) // ( boundary(i) && boundary(j) && i==j )
+ )
+ )
+ {
+ // do nothing, so add entries to interface matrix
+ }
+ else
+ {
+ cell_matrix(i,j) = 0;
+ }
+
- boundary_interface_constraints[cell->level()]
+ empty_constraints
.distribute_local_to_global (cell_matrix,
local_dof_indices,
mg_interface_in[cell->level()]);
template <int dim>
-void LaplaceProblem<dim>::solve ()
+void LaplaceProblem<dim>::solve (bool use_mw)
{
MGTransferPrebuilt<Vector<double> > mg_transfer(hanging_node_constraints, mg_constrained_dofs);
mg_transfer.build_matrices(mg_dof_handler);
mg::Matrix<> mg_matrix(mg_matrices);
mg::Matrix<> mg_interface_up(mg_interface_in);
mg::Matrix<> mg_interface_down(mg_interface_in);
+ //if (use_mw)
+ //mg_interface_down.initialize(mg_interface_out);
Multigrid<Vector<double> > mg(mg_dof_handler,
mg_matrix,
deallog << std::endl;
assemble_system ();
+ assemble_multigrid (false);
+
+ solve (false);
+
assemble_multigrid (true);
+ solve (true);
- solve ();
// output_results (cycle);
}
}
DEAL:: Number of active cells: 4
DEAL::Number of degrees of freedom: 9 L0: 4 L1: 9
DEAL:: Number of degrees of freedom: 9 (by level: 4, 9)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 1.36931
+DEAL:cg::Convergence step 4 value 8.55703e-14
+DEAL:: 4 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 1.36931
DEAL:cg::Convergence step 4 value 8.55703e-14
DEAL:: 4 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 16
DEAL::Number of degrees of freedom: 25 L0: 4 L1: 9 L2: 25
DEAL:: Number of degrees of freedom: 25 (by level: 4, 9, 25)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.843171
+DEAL:cg::Convergence step 6 value 9.78112e-15
+DEAL:: 6 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.843171
DEAL:cg::Convergence step 6 value 9.78112e-15
DEAL:: 6 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 28
DEAL::Number of degrees of freedom: 41 L0: 4 L1: 9 L2: 25 L3: 25
DEAL:: Number of degrees of freedom: 41 (by level: 4, 9, 25, 25)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.737526
+DEAL:cg::Convergence step 8 value 8.81541e-15
+DEAL:: 8 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.737526
DEAL:cg::Convergence step 8 value 8.81541e-15
DEAL:: 8 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 40
DEAL::Number of degrees of freedom: 57 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25
DEAL:: Number of degrees of freedom: 57 (by level: 4, 9, 25, 25, 25)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.730417
+DEAL:cg::Convergence step 8 value 9.36194e-14
+DEAL:: 8 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.730417
DEAL:cg::Convergence step 8 value 9.36194e-14
DEAL:: 8 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 52
DEAL::Number of degrees of freedom: 73 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25 L5: 25
DEAL:: Number of degrees of freedom: 73 (by level: 4, 9, 25, 25, 25, 25)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.729970
+DEAL:cg::Convergence step 8 value 2.81558e-13
+DEAL:: 8 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.729970
DEAL:cg::Convergence step 8 value 2.81558e-13
DEAL:: 8 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 160
DEAL::Number of degrees of freedom: 201 L0: 4 L1: 9 L2: 25 L3: 65 L4: 65 L5: 65 L6: 65
DEAL:: Number of degrees of freedom: 201 (by level: 4, 9, 25, 65, 65, 65, 65)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.492333
+DEAL:cg::Convergence step 9 value 5.07798e-13
+DEAL:: 9 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.492333
DEAL:cg::Convergence step 9 value 5.07798e-13
DEAL:: 9 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 304
DEAL::Number of degrees of freedom: 357 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 81 L7: 153
DEAL:: Number of degrees of freedom: 357 (by level: 4, 9, 25, 81, 81, 81, 81, 153)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.404672
+DEAL:cg::Convergence step 8 value 5.14673e-13
+DEAL:: 8 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.404672
DEAL:cg::Convergence step 8 value 5.14673e-13
DEAL:: 8 CG iterations needed to obtain convergence.
DEAL:: Number of active cells: 676
DEAL::Number of degrees of freedom: 761 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 121 L7: 209 L8: 465
DEAL:: Number of degrees of freedom: 761 (by level: 4, 9, 25, 81, 81, 81, 121, 209, 465)
+DEAL::assemble_multigrid
+DEAL:cg::Starting value 0.404553
+DEAL:cg::Convergence step 9 value 6.34913e-14
+DEAL:: 9 CG iterations needed to obtain convergence.
+DEAL::assemble_multigrid (mesh_worker)
DEAL:cg::Starting value 0.404553
DEAL:cg::Convergence step 9 value 6.34913e-14
DEAL:: 9 CG iterations needed to obtain convergence.
std::vector<std::vector<bool> > interface_dofs
= mg_constrained_dofs.get_refinement_edge_indices ();
- std::vector<std::vector<bool> > boundary_interface_dofs
- = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
- std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
+ ConstraintMatrix empty_constraints;
for (unsigned int level=0; level<triangulation.n_levels(); ++level)
{
boundary_constraints[level].add_lines (interface_dofs[level]);
boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
boundary_constraints[level].close ();
-
- boundary_interface_constraints[level]
- .add_lines (boundary_interface_dofs[level]);
- boundary_interface_constraints[level].close ();
}
typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
// (in the <code>solve()</code>
// function) be able to just pass the
// transpose matrix where necessary.
+ const unsigned int lvl = cell->level();
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ( !(interface_dofs[cell->level()][local_dof_indices[i]]==true &&
- interface_dofs[cell->level()][local_dof_indices[j]]==false))
- cell_matrix(i,j) = 0;
-
- boundary_interface_constraints[cell->level()]
+ if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i)
+ &&
+ !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j)
+ &&
+ (
+ (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j])
+ ) // ( !boundary(i) && !boundary(j) )
+ ||
+ (
+ mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i])
+ &&
+ local_dof_indices[i]==local_dof_indices[j]
+ ) // ( boundary(i) && boundary(j) && i==j )
+ )
+ )
+ {
+ // do nothing, so add entries to interface matrix
+ }
+ else
+ {
+ cell_matrix(i,j) = 0;
+ std::cout << i << " " << j << "\n";
+ }
+
+
+ empty_constraints
.distribute_local_to_global (cell_matrix,
local_dof_indices,
mg_interface_matrices[cell->level()]);
assemble_multigrid ();
solve ();
-// output_results (cycle);
+ output_results (cycle);
}
}
projects / dealii.git / commitdiff