void assemble(MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
- const std::vector<unsigned int>& i2,
- const bool transpose = false);
+ const std::vector<unsigned int>& i2);
+
+ /**
+ * Assemble a single matrix
+ * into a global matrix.
+ */
+ void assemble_transpose(MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2);
/**
* The global matrix being
template <class MATRIX>
inline void
- MGMatrixSimple<MATRIX>::assemble(MATRIX& G,
- const FullMatrix<double>& M,
- const std::vector<unsigned int>& i1,
- const std::vector<unsigned int>& i2,
- bool transpose)
+ MGMatrixSimple<MATRIX>::assemble(
+ MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2)
{
AssertDimension(M.m(), i1.size());
AssertDimension(M.n(), i2.size());
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(j,k)) >= threshold)
- {
- if (transpose)
- G.add(i1[j], i2[k], M(j,k));
- else
- G.add(i1[j], i2[k], M(j,k));
- }
+ G.add(i1[j], i2[k], M(j,k));
+ }
+
+
+ template <class MATRIX>
+ inline void
+ MGMatrixSimple<MATRIX>::assemble_transpose(
+ MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2)
+ {
+ AssertDimension(M.n(), i1.size());
+ AssertDimension(M.m(), i2.size());
+
+ for (unsigned int j=0; j<i1.size(); ++j)
+ for (unsigned int k=0; k<i2.size(); ++k)
+ if (std::fabs(M(k,j)) >= threshold)
+ G.add(i1[j], i2[k], M(k,j));
}
// which is done by
// the coarser cell
assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices);
- assemble((*flux_up)[level1],info1.M2[0].matrix, info1.indices, info2.indices, true);
+ assemble_transpose((*flux_up)[level1],info1.M2[0].matrix, info2.indices, info1.indices);
assemble((*flux_down)[level1], info2.M2[0].matrix, info2.indices, info1.indices);
}
}
// Include files for FiniteElement
// classes and DoFHandler.
#include <fe/fe_q.h>
+#include <fe/fe_dgp.h>
#include <fe/fe_dgq.h>
#include <dofs/dof_tools.h>
#include <multigrid/mg_dof_handler.h>
FullMatrix<double>& matrix_v2u2 = info2.M1[0].matrix;
const unsigned int deg = fe1.get_fe().tensor_degree();
- const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
- const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
- const double penalty = penalty1 + penalty2;
+ double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
+ double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
+ if (info1.cell->has_children() ^ info2.cell->has_children())
+ {
+ Assert (info1.face == info2.face, ExcInternalError());
+ Assert (info1.face->has_children(), ExcInternalError());
+ Assert (info1.cell->has_children(), ExcInternalError());
+ penalty1 *= 2;
+ }
+const double penalty = penalty1 + penalty2;
for (unsigned k=0;k<fe1.n_quadrature_points;++k)
for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
void refine_grid ();
void output_results (const unsigned int cycle) const;
- Triangulation<dim> triangulation;
- const MappingQ1<dim> mapping;
+ Triangulation<dim> triangulation;
+ const MappingQ1<dim> mapping;
const FiniteElement<dim>& fe;
- MGDoFHandler<dim> dof_handler;
+ MGDoFHandler<dim> mg_dof_handler;
+ DoFHandler<dim>& dof_handler;
SparsityPattern sparsity;
SparseMatrix<double> matrix;
Step39<dim>::Step39(const FiniteElement<dim>& fe)
:
fe(fe),
- dof_handler(triangulation)
+ mg_dof_handler(triangulation),
+ dof_handler(mg_dof_handler)
{
GridGenerator::hyper_L(triangulation, -1, 1);
+ triangulation.begin()->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
}
unsigned int n_dofs = dof_handler.n_dofs();
CompressedSparsityPattern c_sparsity(n_dofs);
- const DoFHandler<dim>& dof = dof_handler;
- DoFTools::make_flux_sparsity_pattern(dof, c_sparsity);
+ DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
sparsity.copy_from(c_sparsity);
-
matrix.reinit(sparsity);
solution.reinit(n_dofs);
for (unsigned int level=mg_sparsity.get_minlevel();
level<=mg_sparsity.get_maxlevel();++level)
{
- CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level));
+ CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level));
CompressedSparsityPattern ci_sparsity;
if (level>0)
- ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level));
+ ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level));
- MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level);
+ MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level);
if (level>0)
- MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level);
+ MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level);
mg_sparsity[level].copy_from(c_sparsity);
mg_matrix[level].reinit(mg_sparsity[level]);
MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
info_box.initialize(integrator, fe, mapping);
MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator);
+
+// matrix.print_formatted(std::cout, 2, false, 5, "*");
}
const MatrixIntegrator<dim> local;
MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> >, MatrixIntegrator<dim> >
integrator(local);
- const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+ const unsigned int n_gauss_points = mg_dof_handler.get_fe().tensor_degree()+1;
integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
UpdateFlags update_flags = update_values | update_gradients;
integrator.add_update_flags(update_flags, true, true, true, true);
integrator.initialize(mg_matrix);
- MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
+ integrator.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
+ MeshWorker::IntegrationInfoBox<dim> info_box(mg_dof_handler);
info_box.initialize(integrator, fe, mapping);
- MeshWorker::integration_loop(dof_handler.begin(), dof_handler.end(), info_box, integrator);
+ MeshWorker::integration_loop(mg_dof_handler.begin(), mg_dof_handler.end(), info_box, integrator);
+
+// for (unsigned int l=0;l<triangulation.n_levels();++l)
+// {
+// std::cout << "level matrix " << l << std::endl;
+// mg_matrix[l].print_formatted(std::cout, 2, false, 5, "*");
+// }
+// for (unsigned int l=1;l<triangulation.n_levels();++l)
+// {
+// std::cout << "up matrix " << l << std::endl;
+// mg_matrix_dg_up[l].print_formatted(std::cout, 2, false, 5, "*");
+// std::cout << "down matrix " << l << std::endl;
+// mg_matrix_dg_down[l].print_formatted(std::cout, 2, false, 5, "*");
+// }
}
GrowingVectorMemory<Vector<double> > mem;
MGTransferPrebuilt<Vector<double> > mg_transfer;
- mg_transfer.build_matrices(dof_handler);
+ mg_transfer.build_matrices(mg_dof_handler);
FullMatrix<double> coarse_matrix;
coarse_matrix.copy_from (mg_matrix[0]);
MGCoarseGridHouseholder<double, Vector<double> > mg_coarse;
// Now, we are ready to set up the
// V-cycle operator and the
// multilevel preconditioner.
- Multigrid<Vector<double> > mg(dof_handler, mgmatrix,
+ Multigrid<Vector<double> > mg(mg_dof_handler, mgmatrix,
mg_coarse, mg_transfer,
mg_smoother, mg_smoother);
mg.set_edge_matrices(mgdown, mgup);
PreconditionMG<dim, Vector<double>,
MGTransferPrebuilt<Vector<double> > >
- preconditioner(dof_handler, mg, mg_transfer);
+ preconditioner(mg_dof_handler, mg, mg_transfer);
cg.solve(matrix, solution, right_hand_side, preconditioner);
}
void
Step39<dim>::run(unsigned int n_steps)
{
+ deallog << "Element: " << fe.get_name() << std::endl;
for (unsigned int s=0;s<n_steps;++s)
{
deallog << "Step " << s << std::endl;
<< triangulation.n_levels() << " levels" << std::endl;
setup_system();
- deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, levels";
+ deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, level dofs";
for (unsigned int l=0;l<triangulation.n_levels();++l)
- deallog << ' ' << dof_handler.n_dofs(l);
+ deallog << ' ' << mg_dof_handler.n_dofs(l);
deallog << std::endl;
deallog << "Assemble matrix" << std::endl;
assemble_mg_matrix();
deallog << "Assemble right hand side" << std::endl;
assemble_right_hand_side();
-
+ deallog << "Solve" << std::endl;
solve();
+ deallog << "Error" << std::endl;
error();
output_results(s);
}
int main()
{
- FE_DGQ<2> dgq1(1);
- Step39<2> test1(dgq1);
+ FE_DGQ<2> fe1(1);
+ Step39<2> test1(fe1);
test1.run(7);
}