std::ofstream out ("grid-1.eps");
GridOut grid_out;
grid_out.write_eps (triangulation, out);
-};
+}
// default object, over which the
// triangulation has full control.
triangulation.set_boundary (0);
-};
+}
{
first_grid ();
second_grid ();
-};
+}
std::cout << std::endl;
};
-};
+}
// The following, second function also
std::cout << std::endl;
};
-};
+}
// The following main function just
compute_pi_by_perimeter<2> ();
return 0;
-};
+}
<< std::endl
<< "============================"
<< std::endl;
-};
+}
// ``copy_from'' call.
sparsity_pattern.copy_from (csp);
system_matrix.reinit (sparsity_pattern);
-};
+}
output_table.add_value ("cells", triangulation.n_active_cells());
output_table.add_value ("|u|_1", norm);
output_table.add_value ("error", std::fabs(norm-std::sqrt(3.14159265358/2)));
-};
+}
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
-};
+}
output_table.set_precision("error", 6);
output_table.write_text (std::cout);
std::cout << std::endl;
-};
+}
};
return 0;
-};
+}
beta_function (),
rhs_function (),
boundary_function ()
-{};
+{}
// @sect4{Function: assemble_cell_term}
DGMethod<dim>::~DGMethod ()
{
dof_handler.clear ();
-};
+}
template <int dim>
solution1.reinit (dof_handler.n_dofs());
solution2.reinit (dof_handler.n_dofs());
right_hand_side.reinit (dof_handler.n_dofs());
-};
+}
// @sect4{Function: assemble_system1}
for (unsigned int i=0; i<dofs_per_cell; ++i)
right_hand_side(dofs[i]) += cell_vector(i);
}
-};
+}
// @sect4{Function: assemble_system2}
for (unsigned int i=0; i<dofs_per_cell; ++i)
right_hand_side(dofs[i]) += cell_vector(i);
}
-};
+}
// @sect3{All the rest}
// ready to start the linear solver.
solver.solve (system_matrix, solution, right_hand_side,
preconditioner);
-};
+}
// We refine the grid according to a
data_out.build_patches ();
data_out.write_gnuplot(gnuplot_output);
-};
+}
// The following ``run'' function is
};
return 0;
-};
+}
// straightforward:
template <int dim>
EvaluationBase<dim>::~EvaluationBase ()
- {};
+ {}
EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
{
refinement_cycle = step;
- };
+ }
// @sect4{Point evaluation}
:
evaluation_point (evaluation_point),
results_table (results_table)
- {};
+ {}
// into the table of results:
results_table.add_value ("DoFs", dof_handler.n_dofs());
results_table.add_value ("u(x_0)", point_value);
- };
+ }
:
output_name_base (output_name_base),
output_format (output_format)
- {};
+ {}
// After the description above, the
#endif
data_out.write (out, output_format);
- };
+ }
// example, that much shall be
// sufficient, so we close the
// namespace.
-};
+}
// @sect3{The Laplace solver classes}
Base<dim>::Base (Triangulation<dim> &coarse_grid)
:
triangulation (&coarse_grid)
- {};
+ {}
template <int dim>
Base<dim>::~Base ()
- {};
+ {}
// @sect4{A general solver class}
quadrature (&quadrature),
dof_handler (triangulation),
boundary_values (&boundary_values)
- {};
+ {}
// The destructor is simple, it
Solver<dim>::~Solver ()
{
dof_handler.clear ();
- };
+ }
// The next function is the one
LinearSystem linear_system (dof_handler);
assemble_linear_system (linear_system);
linear_system.solve (solution);
- };
+ }
// As stated above, the
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
{
postprocessor (dof_handler, solution);
- };
+ }
// The ``n_dofs'' function should
Solver<dim>::n_dofs () const
{
return dof_handler.n_dofs();
- };
+ }
// The following function assembles
solution,
linear_system.rhs);
- };
+ }
// The second of this pair of
cell_matrix(i,j));
mutex.release ();
};
- };
+ }
// Now for the functions that
sparsity_pattern.compress();
matrix.reinit (sparsity_pattern);
rhs.reinit (dof_handler.n_dofs());
- };
+ }
cg.solve (matrix, solution, rhs, preconditioner);
hanging_node_constraints.distribute (solution);
- };
+ }
Solver<dim> (triangulation, fe,
quadrature, boundary_values),
rhs_function (&rhs_function)
- {};
+ {}
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs(local_dof_indices[i]) += cell_rhs(i);
};
- };
+ }
// @sect4{Global refinement}
Base<dim> (coarse_grid),
PrimalSolver<dim> (coarse_grid, fe, quadrature,
rhs_function, boundary_values)
- {};
+ {}
RefinementGlobal<dim>::refine_grid ()
{
this->triangulation->refine_global (1);
- };
+ }
// @sect4{Local refinement by the Kelly error indicator}
Base<dim> (coarse_grid),
PrimalSolver<dim> (coarse_grid, fe, quadrature,
rhs_function, boundary_values)
- {};
+ {}
estimated_error_per_cell,
0.3, 0.03);
this->triangulation->execute_coarsening_and_refinement ();
- };
+ }
-};
+}
q += std::sin(10*p(i)+5*p(0)*p(0));
const double exponential = std::exp(q);
return exponential;
-};
+}
t1 = t1*t1;
return -u*(t1+t2+t3);
-};
+}
// Finally end the line in which we
// displayed status reports:
std::cout << std::endl;
-};
+}
// And one blank line after all
// results:
std::cout << std::endl;
-};
+}
};
return 0;
-};
+}
template <int dim>
EvaluationBase<dim>::~EvaluationBase ()
- {};
+ {}
EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
{
refinement_cycle = step;
- };
+ }
// @sect4{The PointValueEvaluation class}
PointValueEvaluation (const Point<dim> &evaluation_point)
:
evaluation_point (evaluation_point)
- {};
+ {}
std::cout << " Point value=" << point_value
<< std::endl;
- };
+ }
// @sect4{The PointXDerivativeEvaluation class}
PointXDerivativeEvaluation (const Point<dim> &evaluation_point)
:
evaluation_point (evaluation_point)
- {};
+ {}
// The more interesting things
point_derivative /= evaluation_point_hits;
std::cout << " Point x-derivative=" << point_derivative
<< std::endl;
- };
+ }
GridOutput (const std::string &output_name_base)
:
output_name_base (output_name_base)
- {};
+ {}
template <int dim>
#endif
GridOut().write_eps (dof_handler.get_tria(), out);
- };
-};
+ }
+}
// @sect3{The Laplace solver classes}
Base<dim>::Base (Triangulation<dim> &coarse_grid)
:
triangulation (&coarse_grid)
- {};
+ {}
template <int dim>
Base<dim>::~Base ()
- {};
+ {}
Base<dim>::set_refinement_cycle (const unsigned int cycle)
{
refinement_cycle = cycle;
- };
+ }
// @sect4{The Laplace Solver class}
face_quadrature (&face_quadrature),
dof_handler (triangulation),
boundary_values (&boundary_values)
- {};
+ {}
template <int dim>
Solver<dim>::~Solver ()
{
dof_handler.clear ();
- };
+ }
template <int dim>
LinearSystem linear_system (dof_handler);
assemble_linear_system (linear_system);
linear_system.solve (solution);
- };
+ }
template <int dim>
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
{
postprocessor (dof_handler, solution);
- };
+ }
template <int dim>
Solver<dim>::n_dofs () const
{
return dof_handler.n_dofs();
- };
+ }
template <int dim>
linear_system.matrix,
solution,
linear_system.rhs);
- };
+ }
template <int dim>
cell_matrix(i,j));
mutex.release ();
};
- };
+ }
template <int dim>
sparsity_pattern.compress();
matrix.reinit (sparsity_pattern);
rhs.reinit (dof_handler.n_dofs());
- };
+ }
cg.solve (matrix, solution, rhs, preconditioner);
hanging_node_constraints.distribute (solution);
- };
+ }
quadrature, face_quadrature,
boundary_values),
rhs_function (&rhs_function)
- {};
+ {}
template <int dim>
PrimalSolver<dim>::solve_problem ()
{
Solver<dim>::solve_problem ();
- };
+ }
PrimalSolver<dim>::n_dofs() const
{
return Solver<dim>::n_dofs();
- };
+ }
template <int dim>
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
{
Solver<dim>::postprocess(postprocessor);
- };
+ }
template <int dim>
#endif
data_out.write (out, DataOut<dim>::gnuplot);
- };
+ }
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs(local_dof_indices[i]) += cell_rhs(i);
};
- };
+ }
// @sect4{The RefinementGlobal and RefinementKelly classes}
PrimalSolver<dim> (coarse_grid, fe, quadrature,
face_quadrature, rhs_function,
boundary_values)
- {};
+ {}
RefinementGlobal<dim>::refine_grid ()
{
this->triangulation->refine_global (1);
- };
+ }
PrimalSolver<dim> (coarse_grid, fe, quadrature,
face_quadrature,
rhs_function, boundary_values)
- {};
+ {}
estimated_error_per_cell,
0.3, 0.03);
this->triangulation->execute_coarsening_and_refinement ();
- };
+ }
face_quadrature,
rhs_function, boundary_values),
weighting_function (&weighting_function)
- {};
+ {}
estimated_error,
0.3, 0.03);
this->triangulation->execute_coarsening_and_refinement ();
- };
+ }
-};
+}
// @sect3{Equation data}
SetUp<Traits,dim>::get_boundary_values () const
{
return boundary_values;
- };
+ }
template <class Traits, int dim>
SetUp<Traits,dim>::get_right_hand_side () const
{
return right_hand_side;
- };
+ }
template <class Traits, int dim>
create_coarse_grid (Triangulation<dim> &coarse_grid) const
{
Traits::create_coarse_grid (coarse_grid);
- };
+ }
// @sect4{The CurvedRidges class}
q += std::sin(10*p(i)+5*p(0)*p(0));
const double exponential = std::exp(q);
return exponential;
- };
+ }
t1 = t1*t1;
return -u*(t1+t2+t3);
- };
+ }
template <int dim>
{
GridGenerator::hyper_cube (coarse_grid, -1, 1);
coarse_grid.refine_global (2);
- };
+ }
// @sect4{The Exercise_2_3 class}
// this example is a grid point,
// we refine once globally:
coarse_grid.refine_global (1);
- };
-};
+ }
+}
// @sect4{Discussion}
//
PointValueEvaluation (const Point<dim> &evaluation_point)
:
evaluation_point (evaluation_point)
- {};
+ {}
// As for doing the main purpose of
// point, so raise an exception
// unconditionally:
AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point));
- };
+ }
// @sect4{The PointValueEvaluation class}
PointXDerivativeEvaluation (const Point<dim> &evaluation_point)
:
evaluation_point (evaluation_point)
- {};
+ {}
// What is interesting is the
// measure of the volume over
// which we have integrated:
rhs.scale (1./total_volume);
- };
+ }
-};
+}
// @sect3{Extending the LaplaceSolver namespace}
quadrature, face_quadrature,
boundary_values),
dual_functional (&dual_functional)
- {};
+ {}
template <int dim>
DualSolver<dim>::solve_problem ()
{
Solver<dim>::solve_problem ();
- };
+ }
DualSolver<dim>::n_dofs() const
{
return Solver<dim>::n_dofs();
- };
+ }
template <int dim>
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
{
Solver<dim>::postprocess(postprocessor);
- };
+ }
assemble_rhs (Vector<double> &rhs) const
{
dual_functional->assemble_rhs (this->dof_handler, rhs);
- };
+ }
// @sect4{The WeightedResidual class}
rhs_values.resize(n_q_points);
dual_weights.resize(n_q_points);
cell_grad_grads.resize(n_q_points);
- };
+ }
dual_weights.resize(n_face_q_points);
cell_grads.resize(n_face_q_points);
neighbor_grads.resize(n_face_q_points);
- };
+ }
DualSolver<dim> (coarse_grid, dual_fe,
quadrature, face_quadrature,
dual_functional)
- {};
+ {}
// The next five functions are
Threads::encapsulate (&WeightedResidual<2>::solve_dual_problem)
.collect_args (this));
thread_manager.wait ();
- };
+ }
template <int dim>
WeightedResidual<dim>::solve_primal_problem ()
{
PrimalSolver<dim>::solve_problem ();
- };
+ }
template <int dim>
void
WeightedResidual<dim>::solve_dual_problem ()
{
DualSolver<dim>::solve_problem ();
- };
+ }
template <int dim>
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
{
PrimalSolver<dim>::postprocess (postprocessor);
- };
+ }
template <int dim>
WeightedResidual<dim>::n_dofs () const
{
return PrimalSolver<dim>::n_dofs();
- };
+ }
error_indicators,
0.8, 0.02);
this->triangulation->execute_coarsening_and_refinement ();
- };
+ }
// Since we want to output both the
#endif
data_out.write (out, DataOut<dim>::gnuplot);
- };
+ }
// @sect3{Estimating errors}
<< std::accumulate (error_indicators.begin(),
error_indicators.end(), 0.)
<< std::endl;
- };
+ }
// @sect4{Estimating on a subset of cells}
if (cell == dual_solver.dof_handler.end())
break;
};
- };
+ }
// @sect4{Computing cell term error contributions}
cell_data.dual_weights[p] *
cell_data.fe_values.JxW (p));
error_indicators(cell_index) += sum;
- };
+ }
// @sect4{Computing edge term error contributions - 1}
// sum up the contributions for
// each cell individually.
face_integrals[cell->face(face_no)] = face_integral;
- };
+ }
// @sect4{Computing edge term error contributions - 2}
// Finally store the value with
// the parent face.
face_integrals[face] = sum;
- };
+ }
-};
+}
// @sect3{A simulation framework}
dual_fe_degree (2),
refinement_criterion (dual_weighted_error_estimator),
max_degrees_of_freedom (20000)
-{};
+{}
std::cout << std::endl;
delete solver;
solver = 0;
-};
+}
};
return 0;
-};
+}
const unsigned int) const
{
return std::pow(p(0), 1./3.);
-};
+}
MinimizationProblem<dim>::MinimizationProblem () :
fe (1),
dof_handler (triangulation)
-{};
+{}
template <int dim>
MinimizationProblem<dim>::~MinimizationProblem ()
{
dof_handler.clear ();
-};
+}
hanging_node_constraints.condense (sparsity_pattern);
sparsity_pattern.compress();
-};
+}
template <int dim>
for (unsigned int k=0; k<n; k+=2)
p += (v*v);
return p;
-};
+}
template <int dim>
newton_update,
newton_residual);
// std::cout << " res=" << newton_residual.l2_norm() << std::endl;
-};
+}
present_solution += newton_update;
hanging_node_constraints.distribute (present_solution);
-};
+}
VectorTools::interpolate (dof_handler,
InitializationValues(),
present_solution);
-};
+}
Vector<double> tmp (dof_handler.n_dofs());
solution_transfer.interpolate (present_solution, tmp);
present_solution = tmp;
-};
+}
};
return energy;
-};
+}
#endif
data_out.write_gnuplot (out);
-};
+}
refine_grid ();
};
-};
+}
int main ()
return 1;
};
return 0;
-};
+}
triangulation.execute_coarsening_and_refinement ();
};
-};
+}
// Up to now, we only have a grid,
// all distributed symmetrically
// around the origin, this shows up
// again in the sparsity pattern.
-};
+}
// 1260 rows and columns, but large
// matrices often have several
// 100,000s).
-};
+}
// of degrees of freedom to the
// sparsity pattern of the matrix.
renumber_dofs (dof_handler);
-};
+}
LaplaceProblem::LaplaceProblem () :
fe (1),
dof_handler (triangulation)
-{};
+{}
// Now, the first thing we've got to
// values:
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
// Now comes the difficult part:
system_matrix,
solution,
system_rhs);
-};
+}
// The following function simply
// job, the solution variable
// contains the nodal values of the
// solution function.
-};
+}
// The last part of a typical finite
// GMV, or some other format):
std::ofstream output ("solution.gpl");
data_out.write_gnuplot (output);
-};
+}
// The following function is the main
assemble_system ();
solve ();
output_results ();
-};
+}
LaplaceProblem laplace_problem;
laplace_problem.run ();
return 0;
-};
+}
return_value += 4*std::pow(p(i), 4);
return return_value;
-};
+}
// The boundary values were to be
const unsigned int) const
{
return p.square();
-};
+}
LaplaceProblem<dim>::LaplaceProblem () :
fe (1),
dof_handler (triangulation)
-{};
+{}
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
system_matrix,
solution,
system_rhs);
-};
+}
// Solving the linear system of
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
-};
+}
"solution-2d.gmv" :
"solution-3d.gmv");
data_out.write_gmv (output);
-};
+}
assemble_system ();
solve ();
output_results ();
-};
+}
laplace_problem_3d.run ();
return 0;
-};
+}
return 20;
else
return 1;
-};
+}
else
values[i] = 1;
};
-};
+}
// This function is as before.
LaplaceProblem<dim>::LaplaceProblem () :
fe (1),
dof_handler (triangulation)
-{};
+{}
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
system_matrix,
solution,
system_rhs);
-};
+}
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
-};
+}
// And then write the data to the
// file.
data_out.write_eps (output);
-};
+}
solve ();
output_results (cycle);
};
-};
+}
*/
return 0;
-};
+}
return 20;
else
return 1;
-};
+}
else
values[i] = 1;
};
-};
+}
// This is mostly the same as before,
LaplaceProblem<dim>::LaplaceProblem () :
fe (2),
dof_handler (triangulation)
-{};
+{}
// Here comes the added destructor of
LaplaceProblem<dim>::~LaplaceProblem ()
{
dof_handler.clear ();
-};
+}
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
system_matrix,
solution,
system_rhs);
-};
+}
// the solutions of the linear
// system just solved.
hanging_node_constraints.distribute (solution);
-};
+}
// Instead of global refinement, we
// perform the actions for which
// the cells are flagged.
triangulation.execute_coarsening_and_refinement ();
-};
+}
// each grid as a postscript file.
GridOut grid_out;
grid_out.write_eps (triangulation, output);
-};
+}
std::ofstream output ("final-solution.eps");
data_out.write_eps (output);
-};
+}
// The main function is unaltered in
// performed as was expected and we
// can return without error.
return 0;
-};
+}
};
return return_value;
-};
+}
};
return return_value;
-};
+}
};
return return_value;
-};
+}
// @sect3{The Laplace solver class}
dof_handler (triangulation),
fe (&fe),
refinement_mode (refinement_mode)
-{};
+{}
LaplaceProblem<dim>::~LaplaceProblem ()
{
dof_handler.clear ();
-};
+}
// The following function sets up the
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
system_matrix,
solution,
system_rhs);
-};
+}
// Solving the system of equations is
preconditioner);
hanging_node_constraints.distribute (solution);
-};
+}
// Now for the function doing grid
break;
};
};
-};
+}
// used:
convergence_table.set_tex_format("cells", "r");
convergence_table.set_tex_format("dofs", "r");
-};
+}
convergence_table.write_tex(table_file);
table_file.close();
}
-};
+}
// @sect3{Main function}
};
return 0;
-};
+}
template <int dim>
RightHandSide<dim>::RightHandSide () :
Function<dim> (dim)
-{};
+{}
// This is the function that returns
values(1) = 1;
else
values(1) = 0;
-};
+}
// resulting code can in principle
// be as fast as if we had
// duplicated the code.
-};
+}
// constructor only needs the
// parameter to generate a copy of
// the finite element from this.
-{};
+{}
ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
-};
+}
// Setting up the system of equations
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
// The big changes in this program
system_matrix,
solution,
system_rhs);
-};
+}
preconditioner);
hanging_node_constraints.distribute (solution);
-};
+}
0.3, 0.03);
triangulation.execute_coarsening_and_refinement ();
-};
+}
// The output happens mostly as has
data_out.add_data_vector (solution, solution_names);
data_out.build_patches ();
data_out.write_gmv (output);
-};
+}
solve ();
output_results (cycle);
};
-};
+}
// The main function is again exactly
};
return 0;
-};
+}
value[i] = 1+0.8*std::sin(8*M_PI*p[0]);
return value;
-};
+}
for (unsigned int i=0; i<points.size(); ++i)
values[i] = AdvectionField<dim>::value (points[i]);
-};
+}
return ( (p-center_point).square() < diameter*diameter ?
.1/std::pow(diameter,dim) :
0);
-};
+}
for (unsigned int i=0; i<points.size(); ++i)
values[i] = RightHandSide<dim>::value (points[i], component);
-};
+}
const double sine_term = std::sin(16*M_PI*std::sqrt(p.square()));
const double weight = std::exp(-5*p.square()) / std::exp(-5.);
return sine_term * weight;
-};
+}
for (unsigned int i=0; i<points.size(); ++i)
values[i] = BoundaryValues<dim>::value (points[i], component);
-};
+}
AdvectionProblem<dim>::AdvectionProblem () :
dof_handler (triangulation),
fe(1)
-{};
+{}
AdvectionProblem<dim>::~AdvectionProblem ()
{
dof_handler.clear ();
-};
+}
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-};
+}
// due to the fact that we have
// included them into the weak
// formulation of the problem.
-};
+}
// therefore also operate on
// the same ``lock''.
};
-};
+}
preconditioner);
hanging_node_constraints.distribute (solution);
-};
+}
// The following function refines the
0.5, 0.03);
triangulation.execute_coarsening_and_refinement ();
-};
+}
GridOut grid_out;
grid_out.write_eps (triangulation, output);
-};
+}
// ... as is the main loop (setup -
std::ofstream output ("final-solution.gmv");
data_out.write_gmv (output);
-};
+}
// eliminating the need to write
// code included in conditional
// preprocessor sections.
-};
+}
// Following now the function that
1+1.0*dim/2) *
std::sqrt(gradient.square()));
};
-};
+}
// @sect3{Main function}
};
return 0;
-};
+}