// Delete values on refinement edge,
// since restriction will add them again.
for (unsigned int i=0; i<dofs_per_face; ++i)
- defect[level](level_face_indices[i]) = 0.;
+ defect[level](level_face_indices[i])
+ = 0.;
};
};
};
const Vector<double>& v,
const char* name) const
{
- Vector<double> out_vector;
-
const DoFHandler<dim>* dof = mg_dof_handler;
+ Vector<double> out_vector(dof->n_dofs());
+
+ out_vector = -10000;
+
const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
vector<unsigned int> global_dof_indices (dofs_per_cell);
vector<unsigned int> level_dof_indices (dofs_per_cell);
- DoFHandler<dim>::cell_iterator
- global_cell = dof->begin(level);
- MGDoFHandler<dim>::cell_iterator
- level_cell = mg_dof_handler->begin(level);
+ DoFHandler<dim>::active_cell_iterator
+ global_cell = dof->begin_active(level);
+ MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler->begin_active(level);
const MGDoFHandler<dim>::cell_iterator
endc = mg_dof_handler->end(level);
const MGCoarseGridSolver &coarse_grid_solver)
{
level_mgstep (maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
- abort ();
+// abort ();
}
if (level == minlevel)
{
coarse_grid_solver(level, solution[level], defect[level]);
+ sprintf(name, "MG%d-solution",level);
+ print_vector(level, solution[level], name);
return;
}
solution[level] += t;
- // smoothing (modify solution again)
-//TODO: what happens here? smooth overwrites the solution[level],
-//TODO: so the previous two statements should have no effect. No?
+ // post-smoothing
+
post_smooth.smooth(level, solution[level], defect[level]);
sprintf(name, "MG%d-post",level);
print_vector(level, solution[level], name);
+
+ delete[] name;
}
cell->set_refine_flag();
tr.execute_coarsening_and_refinement();
- tr.refine_global(1);
+ tr.refine_global(2);
dof.distribute_dofs(*fe);
const unsigned int size = dof.n_dofs();
deallog << "DoFs " << size << endl;
deallog << "Levels: " << tr.n_levels() << endl;
- for (unsigned int step=1;step < 5; ++step)
+ for (unsigned int step=14;step < 15; ++step)
{
deallog << "smoothing-steps" << step << endl;
SparsityPattern structure(size, dof.max_couplings_between_dofs());
hanging_nodes.condense(A);
hanging_nodes.condense(f);
+ if (true)
+ {
+ ofstream out_file("MGf");
+ DataOut<2> out;
+ out.attach_dof_handler(dof);
+ out.add_data_vector(f, "v");
+ out.build_patches(5);
+ out.write_gnuplot(out_file);
+ }
+
Vector<double> u;
u.reinit(f);
PrimitiveVectorMemory<> mem;
}
equation.build_mgmatrix(mgA, mgdof, quadrature);
- SolverControl cgcontrol(20,0., false, false);
+ SolverControl cgcontrol(2000,1.e-14, false, false);
PrimitiveVectorMemory<> cgmem;
SolverCG<> cgsolver(cgcontrol, cgmem);
PreconditionIdentity cgprec;
hanging_nodes.distribute(u);
DataOut<2> out;
- ofstream ofile("out.gnuplot");
+ char* name = new char[100];
+
+ sprintf(name, "MG-Q%d-%d", degree, step);
+
+ ofstream ofile(name);
out.attach_dof_handler(dof);
out.add_data_vector(u,"u");
- out.build_patches(4);
+ out.add_data_vector(f,"f");
+ out.build_patches(5);
out.write_gnuplot(ofile);
+
+ delete[] name;
}
deallog.pop();
}
RHSFunction<dim>::value (const Point<dim>&p,
const unsigned int) const
{
-// return 1.;
+ return 1.;
return p(0)*p(0)+p(1)*p(1);
return (2.1)*(sin(p(0))* sin(p(1)));
}
-
-
-void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector<double>& v, unsigned int level,
- ostream &out)
-{
- MGDoFHandler<2>::cell_iterator cell;
- MGDoFHandler<2>::cell_iterator endc = dofs.end(level);
-
- Vector<double> values(dofs.get_fe().dofs_per_cell);
-
- unsigned int cell_index=0;
- for (cell=dofs.begin(level); cell!=endc; ++cell, ++cell_index)
- {
- cell->get_mg_dof_values(v, values);
-
- out << cell->vertex(0) << " " << values(0) << endl;
- out << cell->vertex(1) << " " << values(1) << endl;
- out << endl;
- out << cell->vertex(3) << " " << values(3) << endl;
- out << cell->vertex(2) << " " << values(2) << endl;
- out << endl;
- out << endl;
- }
- out << endl;
-}