Vector<double> &dst,
const Vector<double> &src,
const Vector<double> &rhs) = 0;
-
+
+ /**
+ * Print a level vector using
+ * @ref{DoFHandler}.
+ */
+ virtual void print_vector (const unsigned int level,
+ const Vector<double>& v,
+ const char* name) const = 0;
+
private:
/**
* Auxiliary vector.
const Vector<double> &src,
const Vector<double> &rhs);
+ /**
+ * Print a level vector using
+ * @ref{DoFHandler}.
+ */
+ virtual void print_vector (const unsigned int level,
+ const Vector<double>& v,
+ const char* name) const;
+
/**
* Read-only access to level matrices.
*/
#include <dofs/dof_constraints.h>
+#include <numerics/data_out.h>
#include <multigrid/multigrid.h>
#include <algorithm>
#include <fstream>
// traverse all cells and copy the
// data appropriately to the output
// vector
+
+ // Is the level monotonuosly increasing?
+
for (; level_cell != endc; ++level_cell, ++global_cell)
{
const unsigned int level = level_cell->level();
}
+
+template <int dim>
+void
+Multigrid<dim>::print_vector (const unsigned int level,
+ const Vector<double>& v,
+ const char* name) const
+{
+ Vector<double> out_vector;
+
+ const DoFHandler<dim>* dof = mg_dof_handler;
+
+ 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);
+ const MGDoFHandler<dim>::cell_iterator
+ endc = mg_dof_handler->end(level);
+
+ // traverse all cells and copy the
+ // data appropriately to the output
+ // vector
+ for (; level_cell != endc; ++level_cell, ++global_cell)
+ {
+ global_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)
+ out_vector(global_dof_indices[i])
+ = v(level_dof_indices[i]);
+ }
+
+ ofstream out_file(name);
+ DataOut<dim> out;
+ out.attach_dof_handler(*dof);
+ out.add_data_vector(out_vector, "v");
+ out.build_patches(5);
+ out.write_gnuplot(out_file);
+}
+
#endif
+
#include <cmath>
-//TODO: this function is only for debugging purposes and should be removed sometimes
-template<class VECTOR>
-static
-void print_vector(ostream& s, const VECTOR& v, const char* text)
-{
- const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
- unsigned int k=0;
-
- s << endl << "splot '-' title '" << text << "'" << endl;
-
- // write the vector entries in a
- // kind of square
- for (unsigned int i=0;i<n;++i)
- {
- for (unsigned int j=0;j<n;++j)
- s << '\n' << v(k++);
- s << endl;
- }
- s << "e\npause -1" << endl;
-}
-
-
-//////////////////////////////////////////////////////////////////////
-
-
MGBase::~MGBase ()
-{};
+{}
MGBase::MGBase(const MGTransferBase &transfer,
{
Assert(minlevel <= maxlevel,
ExcSwitchedLevels(minlevel, maxlevel));
-};
+}
void
const MGSmootherBase &post_smooth,
const MGCoarseGridSolver &coarse_grid_solver)
{
-// static int k=0;
-// cout << "set title 'cycle " << ++k << "'\n";
-
- level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-};
+ level_mgstep (maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
+ abort ();
+}
void
const MGSmootherBase &post_smooth,
const MGCoarseGridSolver &coarse_grid_solver)
{
+ char *name = new char[100];
+
+ sprintf(name, "MG%d-defect",level);
+ print_vector(level, defect[level], name);
+
solution[level] = 0.;
if (level == minlevel)
// smoothing of the residual by modifying s
pre_smooth.smooth(level, solution[level], defect[level]);
// t = d-As
+
+ sprintf(name, "MG%d-pre",level);
+ print_vector(level, solution[level], name);
+
t.reinit(solution[level].size());
level_vmult(level, t, solution[level], defect[level]);
-// print_vector(cout,t,"T");
// make t rhs of lower level
//TODO: this function adds the restricted t to defect[level-1].
//TODO: why don't we have to clear it before?
transfer->restrict_and_add (level, defect[level-1], t);
-// print_vector(cout,defect[level-1],"Dl-1");
// do recursion
level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
t.reinit(solution[level].size());
// do coarse grid correction
-// print_vector(cout,solution[level-1],"Sl-1");
+
transfer->prolongate(level, t, solution[level-1]);
-// print_vector(cout,t,"T");
+
+ sprintf(name, "MG%d-cgc",level);
+ print_vector(level, t, name);
+
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_smooth.smooth(level, solution[level], defect[level]);
-};
+
+ sprintf(name, "MG%d-post",level);
+ print_vector(level, solution[level], name);
+}
//////////////////////////////////////////////////////////////////////
@$(CXX) $(flags) -c $< -o $@
-all: grid_test.check dof_test.check fe_tables.check gradients.check \
+all: grid_generator.check grid_test.check dof_test.check \
+ fe_tables.check gradients.check \
constraints.check mg.check mglocal.check wave-test-3.check second_derivatives.check
exe: $(all:.check=.testcase) benchmark
@$(CXX) $(flags) -o $@ $^
+############################################################
+
+grid_generator-cc-files = grid_generator.cc
+
+ifeq ($(debug-mode),on)
+grid_generator-o-files = $(grid_generator-cc-files:.cc=.go)
+else
+grid_generator-o-files = $(grid_generator-cc-files:.cc=.o)
+endif
+
+grid_generator.testcase: $(grid_generator-o-files) $(libdeal3d) $(libraries)
+ @echo =====linking======= $<
+ @$(CXX) $(flags) -o $@ $^
+
+
############################################################
second_derivatives-cc-files = second_derivatives.cc
#include <dofs/dof_accessor.h>
#include <multigrid/mg_dof_accessor.h>
#include <grid/grid_generator.h>
-#include <numerics/data_io.h>
+#include <numerics/data_out.h>
#include <fe/fe_lib.lagrange.h>
#include <fe/fe_values.h>
#include <multigrid/multigrid.h>
int main()
{
ofstream logfile("mglocal.output");
- logfile.setf(ios::fixed);
- logfile.precision (3);
+// logfile.setf(ios::fixed);
+// logfile.precision (3);
deallog.attach(logfile);
- deallog.depth_console(0);
+// deallog.depth_console(0);
Helmholtz equation;
RHSFunction<2> rhs;
FEQ2<2> fe2;
FEQ3<2> fe3;
FEQ4<2> fe4;
- for (unsigned int degree=1;degree<=1;degree++)
+ for (unsigned int degree=1;degree<=4;degree++)
{
Triangulation<2> tr;
MGDoFHandler<2> mgdof(tr);
cell->set_refine_flag();
tr.execute_coarsening_and_refinement();
- tr.refine_global(2);
+ tr.refine_global(1);
dof.distribute_dofs(*fe);
const unsigned int size = dof.n_dofs();
deallog << "DoFs " << size << endl;
Vector<double> u;
u.reinit(f);
PrimitiveVectorMemory<> mem;
- SolverControl control(20, 1.e-12, true);
- SolverRichardson<> solver(control, mem);
+ SolverControl control(20, 1.e-8, true);
+ SolverCG<> solver(control, mem);
MGLevelObject<SparsityPattern> mgstruct(0, tr.n_levels()-1);
MGLevelObject<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
solver.solve(A, u, f, mgprecondition);
hanging_nodes.distribute(u);
+
+ DataOut<2> out;
+ ofstream ofile("out.gnuplot");
+ out.attach_dof_handler(dof);
+ out.add_data_vector(u,"u");
+ out.build_patches(4);
+ out.write_gnuplot(ofile);
}
deallog.pop();
}