// cos(2*3.1415926536*p(0)));
return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
case 2:
- return ((1-3.1415926536*3.1415926536) *
- cos(3.1415926536*p(0)) *
- cos(3.1415926536*p(1)));
+// return ((1-3.1415926536*3.1415926536) *
+// cos(3.1415926536*p(0)) *
+// cos(3.1415926536*p(1)));
+ return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1)
+ - 3./2.*(p(0)*p(0)+p(1)*p(1))
+ - 6*(p(0)+p(1))
+ + 6);
default:
return 0;
};
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) /*+
+ fe_values.shape_grad(j,point) +
fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)*/) *
+ fe_values.shape_value(j,point)) *
fe_values.JxW(point);
rhs[0](i) += fe_values.shape_value(i,point) *
right_hand_side(fe_values.quadrature_point(point)) *
// tria.create_hyper_ball(Point<2>(2,3),4);
// tria.set_boundary (&boundary);
- tria.refine_global (1);
+/* tria.refine_global (1);
tria.begin_active()->set_refine_flag();
tria.execute_refinement ();
+ tria.refine_global (3);
+*/
+
+ cout << "Making grid..." << endl;
+ const unsigned int dim=2;
+ tria.refine_global (1);
+
+ Triangulation<dim>::active_cell_iterator cell, endc;
+ for (int i=0; i<8; ++i)
+ {
+ int n_levels = tria.n_levels();
+ cell = tria.begin_active();
+ endc = tria.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ double r = rand()*1.0/RAND_MAX,
+ weight = 1.*
+ (cell->level()*cell->level()) /
+ (n_levels*n_levels);
+
+ if (r <= 0.5*weight)
+ cell->set_refine_flag ();
+ };
+
+ tria.execute_refinement ();
+ };
+ tria.refine_global (1);
+
+ cout << "Distributing dofs..." << endl;
dof.distribute_dofs (fe);
+
+ cout << "Assembling matrices..." << endl;
problem.assemble (equation, quadrature, fe);
-
-// problem.solve ();
+ cout << "Solving..." << endl;
+ problem.solve ();
+ cout << "Printing..." << endl;
DataOut<2> out;
ofstream gnuplot("gnuplot.out.5");
problem.fill_data (out);
out.write_gnuplot (gnuplot);
+ gnuplot.close ();
return 0;
};
#include <lac/dsmatrix.h>
#include <base/exceptions.h>
+#include <grid/dof_constraints.h>
+
// forward declaration
template <int dim> class Triangulation;
*/
dVector solution;
+ /**
+ * List of constraints introduced by
+ * hanging nodes.
+ */
+ ConstraintMatrix constraints;
+
friend class Assembler<dim>;
};
void ConstraintMatrix::condense (dVector &vec) const {
+ Assert (sorted == true, ExcMatrixNotClosed());
+
vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
for (unsigned int row=0; row<(unsigned int)vec.n(); ++row)
if (row == (*next_constraint).line)
+void ConstraintMatrix::distribute (const dVector &condensed,
+ dVector &uncondensed) const {
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert ((unsigned int)condensed.n()+n_constraints() == (unsigned int)uncondensed.n(),
+ ExcWrongDimension());
+
+ // store for each line of the new vector
+ // its old line number before
+ // distribution. If the shift is
+ // -1, this line was condensed away
+ vector<int> old_line;
+
+ old_line.reserve (uncondensed.n());
+
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ unsigned int shift = 0;
+ unsigned int n_rows = (unsigned int)uncondensed.n();
+ for (unsigned int row=0; row!=n_rows; ++row)
+ if (row == (*next_constraint).line)
+ {
+ // this line is constrained
+ old_line.push_back (-1);
+ // note that #lines# is ordered
+ ++next_constraint;
+ ++shift;
+ }
+ else
+ old_line.push_back (row-shift);
+
+
+ next_constraint = lines.begin();
+ for (unsigned int line=0; line<(unsigned int)uncondensed.n(); ++line)
+ if (old_line[line] != -1)
+ // line was not condensed away
+ uncondensed(line) = condensed(old_line[line]);
+ else
+ {
+ // line was condensed away, create it newly
+ // first set it to zero
+ uncondensed(line) = 0.;
+ // then add the different contributions
+ for (unsigned int i=0; i<(*next_constraint).entries.size(); ++i)
+ uncondensed(line) += (condensed(old_line[(*next_constraint).entries[i].first]) *
+ (*next_constraint).entries[i].second);
+ ++next_constraint;
+ };
+};
+
+
+
+void ConstraintMatrix::distribute (dVector &vec) const {
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ for (; next_constraint != lines.end(); ++next_constraint)
+ {
+ // make entry in line next_constraint.line
+ // first set it to zero
+ vec((*next_constraint).line) = 0.;
+ // then add the different contributions
+ for (unsigned int i=0; i<(*next_constraint).entries.size(); ++i)
+ vec((*next_constraint).line) += (vec((*next_constraint).entries[i].first) *
+ (*next_constraint).entries[i].second);
+ };
+};
+
+
+
+
unsigned int ConstraintMatrix::n_constraints () const {
return lines.size();
};
system_sparsity.reinit (dof_handler->n_dofs(),
dof_handler->n_dofs(),
dof_handler->max_couplings_between_dofs());
- solution.reinit (dof_handler->n_dofs());
for (unsigned int i=0; i<right_hand_sides.size(); ++i)
right_hand_sides[i]->reinit (dof_handler->n_dofs());
// make up sparsity pattern and
// compress with constraints
- ConstraintMatrix constraints;
-
+ constraints.clear ();
dof_handler->make_constraint_matrix (constraints);
dof_handler->make_sparsity_pattern (system_sparsity);
constraints.condense (system_sparsity);
// reinite system matrix
system_matrix.reinit (system_sparsity);
+ solution.reinit (dof_handler->n_dofs());
// create assembler
AssemblerData<dim> data (*dof_handler,
true, true, //assemble matrix and rhs
- right_hand_sides.size(), //one rhs
+ right_hand_sides.size(),
*this,
quadrature,
fe);
// condense system matrix in-place
constraints.condense (system_matrix);
+ // condense right hand sides in-place
+ for (unsigned int i=0; i<right_hand_sides.size(); ++i)
+ constraints.condense (*right_hand_sides[i]);
};
template <int dim>
void ProblemBase<dim>::solve () {
int max_iter = 4000;
- double tolerance = 1.e-14;
+ double tolerance = 1.e-16;
Control control1(max_iter,tolerance);
PrimitiveVectorMemory<dVector> memory(right_hand_sides[0]->n());
CG<dSMatrix,dVector> cg(control1,memory);
+ // solve
cg (system_matrix, solution, *right_hand_sides[0]);
+ // distribute solution
+ constraints.distribute (solution);
};
// cos(2*3.1415926536*p(0)));
return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
case 2:
- return ((1-3.1415926536*3.1415926536) *
- cos(3.1415926536*p(0)) *
- cos(3.1415926536*p(1)));
+// return ((1-3.1415926536*3.1415926536) *
+// cos(3.1415926536*p(0)) *
+// cos(3.1415926536*p(1)));
+ return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1)
+ - 3./2.*(p(0)*p(0)+p(1)*p(1))
+ - 6*(p(0)+p(1))
+ + 6);
default:
return 0;
};
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) /*+
+ fe_values.shape_grad(j,point) +
fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)*/) *
+ fe_values.shape_value(j,point)) *
fe_values.JxW(point);
rhs[0](i) += fe_values.shape_value(i,point) *
right_hand_side(fe_values.quadrature_point(point)) *
// tria.create_hyper_ball(Point<2>(2,3),4);
// tria.set_boundary (&boundary);
- tria.refine_global (1);
+/* tria.refine_global (1);
tria.begin_active()->set_refine_flag();
tria.execute_refinement ();
+ tria.refine_global (3);
+*/
+
+ cout << "Making grid..." << endl;
+ const unsigned int dim=2;
+ tria.refine_global (1);
+
+ Triangulation<dim>::active_cell_iterator cell, endc;
+ for (int i=0; i<8; ++i)
+ {
+ int n_levels = tria.n_levels();
+ cell = tria.begin_active();
+ endc = tria.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ double r = rand()*1.0/RAND_MAX,
+ weight = 1.*
+ (cell->level()*cell->level()) /
+ (n_levels*n_levels);
+
+ if (r <= 0.5*weight)
+ cell->set_refine_flag ();
+ };
+
+ tria.execute_refinement ();
+ };
+ tria.refine_global (1);
+
+ cout << "Distributing dofs..." << endl;
dof.distribute_dofs (fe);
+
+ cout << "Assembling matrices..." << endl;
problem.assemble (equation, quadrature, fe);
-
-// problem.solve ();
+ cout << "Solving..." << endl;
+ problem.solve ();
+ cout << "Printing..." << endl;
DataOut<2> out;
ofstream gnuplot("gnuplot.out.5");
problem.fill_data (out);
out.write_gnuplot (gnuplot);
+ gnuplot.close ();
return 0;
};