//---------------------------- step-14.cc ---------------------------
// $Id$
-// Version: $Name$
+// Version: $Name$
//
// Copyright (C) 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
namespace Evaluation
{
template <int dim>
- class EvaluationBase
+ class EvaluationBase
{
public:
virtual ~EvaluationBase ();
void set_refinement_cycle (const unsigned int refinement_cycle);
-
+
virtual void operator () (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const = 0;
protected:
template <int dim>
EvaluationBase<dim>::~EvaluationBase ()
{}
-
-
+
+
template <int dim>
void
EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
{
public:
PointValueEvaluation (const Point<dim> &evaluation_point);
-
+
virtual void operator () (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
-
+
DeclException1 (ExcEvaluationPointNotFound,
Point<dim>,
<< "The evaluation point " << arg1
:
evaluation_point (evaluation_point)
{}
-
+
template <int dim>
void
PointValueEvaluation<dim>::
operator () (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution) const
+ const Vector<double> &solution) const
{
double point_value = 1e20;
{
public:
PointXDerivativeEvaluation (const Point<dim> &evaluation_point);
-
+
virtual void operator () (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
-
+
DeclException1 (ExcEvaluationPointNotFound,
Point<dim>,
<< "The evaluation point " << arg1
:
evaluation_point (evaluation_point)
{}
-
+
template <int dim>
void
PointXDerivativeEvaluation<dim>::
operator () (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution) const
+ const Vector<double> &solution) const
{
double point_derivative = 0;
update_gradients | update_q_points);
std::vector<Tensor<1,dim> >
solution_gradients (vertex_quadrature.n_quadrature_points);
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
}
-
+
template <int dim>
class GridOutput : public EvaluationBase<dim>
{
public:
GridOutput (const std::string &output_name_base);
-
+
virtual void operator () (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
private:
:
output_name_base (output_name_base)
{}
-
+
template <int dim>
void
<< this->refinement_cycle
<< ".eps"
<< std::ends;
-
+
GridOut().write_eps (dof_handler.get_tria(),
deallog.get_file_stream());
}
}
-
+
namespace LaplaceSolver
{
template <int dim> class WeightedResidual;
-
-
+
+
template <int dim>
class Base
virtual void set_refinement_cycle (const unsigned int cycle);
virtual void output_solution () const = 0;
-
+
protected:
const SmartPointer<Triangulation<dim> > triangulation;
template <int dim>
- Base<dim>::~Base ()
+ Base<dim>::~Base ()
{}
{
refinement_cycle = cycle;
}
-
+
template <int dim>
Solver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const Quadrature<dim-1> &face_quadrature,
const Function<dim> &boundary_values);
virtual
~Solver ();
virtual
unsigned int
n_dofs () const;
-
+
protected:
const SmartPointer<const FiniteElement<dim> > fe;
const SmartPointer<const Quadrature<dim> > quadrature;
- const SmartPointer<const Quadrature<dim-1> > face_quadrature;
+ const SmartPointer<const Quadrature<dim-1> > face_quadrature;
DoFHandler<dim> dof_handler;
Vector<double> solution;
const SmartPointer<const Function<dim> > boundary_values;
virtual void assemble_rhs (Vector<double> &rhs) const = 0;
-
+
private:
struct LinearSystem
{
LinearSystem (const DoFHandler<dim> &dof_handler);
void solve (Vector<double> &solution) const;
-
+
ConstraintMatrix hanging_node_constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> matrix;
Base<dim> (triangulation),
fe (&fe),
quadrature (&quadrature),
- face_quadrature (&face_quadrature),
+ face_quadrature (&face_quadrature),
dof_handler (triangulation),
boundary_values (&boundary_values)
{}
template <int dim>
- Solver<dim>::~Solver ()
+ Solver<dim>::~Solver ()
{
dof_handler.clear ();
}
{
return dof_handler.n_dofs();
}
-
+
template <int dim>
void
const unsigned int n_threads = multithread_info.n_default_threads;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
- thread_ranges
+ thread_ranges
= Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
dof_handler.end (),
n_threads);
0,
*boundary_values,
boundary_value_map);
-
+
threads.join_all ();
linear_system.hanging_node_constraints.condense (linear_system.matrix);
const typename DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const
{
- FEValues<dim> fe_values (*fe, *quadrature,
+ FEValues<dim> fe_values (*fe, *quadrature,
update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = fe->dofs_per_cell;
void (*mhnc_p) (const DoFHandler<dim> &,
ConstraintMatrix &)
= &DoFTools::make_hanging_node_constraints;
-
+
Threads::Thread<>
mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints);
virtual
void solve_problem ();
-
+
virtual
unsigned int n_dofs () const;
-
+
virtual
void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
virtual
void output_solution () const;
-
+
protected:
const SmartPointer<const Function<dim> > rhs_function;
virtual void assemble_rhs (Vector<double> &rhs) const;
<< std::ends;
data_out.write (deallog.get_file_stream(), DataOut<dim>::gnuplot);
}
-
+
template <int dim>
void
PrimalSolver<dim>::
- assemble_rhs (Vector<double> &rhs) const
+ assemble_rhs (Vector<double> &rhs) const
{
- FEValues<dim> fe_values (*this->fe, *this->quadrature,
+ FEValues<dim> fe_values (*this->fe, *this->quadrature,
update_values | update_q_points |
update_JxW_values);
rhs_function->value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i,q_point) *
for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
estimated_error(cell_index)
*= weighting_function->value (cell->center());
-
+
GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
estimated_error,
0.3, 0.03);
namespace Data
{
-
+
template <int dim>
struct SetUpBase : public Subscriptor
{
virtual
const Function<dim> & get_right_hand_side () const;
-
+
virtual
void create_coarse_grid (Triangulation<dim> &coarse_grid) const;
template <class Traits, int dim>
const Function<dim> &
- SetUp<Traits,dim>::get_boundary_values () const
+ SetUp<Traits,dim>::get_boundary_values () const
{
return boundary_values;
}
template <class Traits, int dim>
const Function<dim> &
- SetUp<Traits,dim>::get_right_hand_side () const
+ SetUp<Traits,dim>::get_right_hand_side () const
{
return right_hand_side;
}
template <class Traits, int dim>
void
SetUp<Traits,dim>::
- create_coarse_grid (Triangulation<dim> &coarse_grid) const
+ create_coarse_grid (Triangulation<dim> &coarse_grid) const
{
Traits::create_coarse_grid (coarse_grid);
}
-
+
template <int dim>
{
public:
BoundaryValues () : Function<dim> () {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component) const;
};
{
public:
RightHandSide () : Function<dim> () {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component) const;
};
void
create_coarse_grid (Triangulation<dim> &coarse_grid);
};
-
-
+
+
template <int dim>
double
CurvedRidges<dim>::BoundaryValues::
100*std::sin(10*p(i)+5*p(0)*p(0));
};
t1 = t1*t1;
-
+
return -u*(t1+t2+t3);
}
GridGenerator::hyper_cube (coarse_grid, -1, 1);
coarse_grid.refine_global (2);
}
-
-
+
+
template <int dim>
struct Exercise_2_3
{
public:
RightHandSide () : ConstantFunction<dim> (1.) {}
};
-
+
static
void
create_coarse_grid (Triangulation<dim> &coarse_grid);
Point<2> (0., -1.),
Point<2> (+1./2, -1.),
Point<2> (+1, -1.),
-
+
Point<2> (-1., -1./2.),
Point<2> (-1./2, -1./2.),
Point<2> (0., -1./2.),
Point<2> (+1./2, -1./2.),
Point<2> (+1, -1./2.),
-
+
Point<2> (-1., 0.),
Point<2> (-1./2, 0.),
Point<2> (+1./2, 0.),
Point<2> (+1, 0.),
-
+
Point<2> (-1., 1./2.),
Point<2> (-1./2, 1./2.),
Point<2> (0., 1./2.),
Point<2> (+1./2, 1./2.),
Point<2> (+1, 1./2.),
-
+
Point<2> (-1., 1.),
Point<2> (-1./2, 1.),
- Point<2> (0., 1.),
+ Point<2> (0., 1.),
Point<2> (+1./2, 1.),
Point<2> (+1, 1.) };
const unsigned int
n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
- for (unsigned int i=0; i<n_cells; ++i)
+ for (unsigned int i=0; i<n_cells; ++i)
{
for (unsigned int j=0;
j<GeometryInfo<dim>::vertices_per_cell;
coarse_grid.create_triangulation (vertices,
cells,
SubCellData());
-
+
coarse_grid.refine_global (1);
}
}
namespace DualFunctional
{
-
+
template <int dim>
class DualFunctionalBase : public Subscriptor
{
};
-
+
template <int dim>
class PointValueEvaluation : public DualFunctionalBase<dim>
{
void
assemble_rhs (const DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const;
-
+
DeclException1 (ExcEvaluationPointNotFound,
Point<dim>,
<< "The evaluation point " << arg1
:
evaluation_point (evaluation_point)
{}
-
+
template <int dim>
void
}
-
+
template <int dim>
class PointXDerivativeEvaluation : public DualFunctionalBase<dim>
{
void
assemble_rhs (const DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const;
-
+
DeclException1 (ExcEvaluationPointNotFound,
Point<dim>,
<< "The evaluation point " << arg1
:
evaluation_point (evaluation_point)
{}
-
+
template <int dim>
void
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
double total_volume = 0;
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
{
fe_values.reinit (cell);
cell_rhs = 0;
-
+
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs.scale (1./total_volume);
}
-
+
}
virtual
void
solve_problem ();
-
+
virtual
unsigned int
n_dofs () const;
{
Solver<dim>::postprocess(postprocessor);
}
-
+
template <int dim>
void
DualSolver<dim>::
- assemble_rhs (Vector<double> &rhs) const
+ assemble_rhs (Vector<double> &rhs) const
{
dual_functional->assemble_rhs (this->dof_handler, rhs);
}
virtual
void
postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
-
+
virtual
unsigned int
n_dofs () const;
private:
void solve_primal_problem ();
void solve_dual_problem ();
-
+
typedef
typename DoFHandler<dim>::active_cell_iterator
active_cell_iterator;
const SmartPointer<const Function<dim> > right_hand_side;
std::vector<double> cell_residual;
- std::vector<double> rhs_values;
- std::vector<double> dual_weights;
+ std::vector<double> rhs_values;
+ std::vector<double> dual_weights;
typename std::vector<Tensor<2,dim> > cell_grad_grads;
CellData (const FiniteElement<dim> &fe,
const Quadrature<dim> &quadrature,
FESubfaceValues<dim> fe_subface_values_cell;
std::vector<double> jump_residual;
- std::vector<double> dual_weights;
+ std::vector<double> dual_weights;
typename std::vector<Tensor<1,dim> > cell_grads;
typename std::vector<Tensor<1,dim> > neighbor_grads;
FaceData (const FiniteElement<dim> &fe,
const Quadrature<dim-1> &face_quadrature);
};
-
+
void estimate_error (Vector<float> &error_indicators) const;
update_q_points |
update_JxW_values),
right_hand_side (&right_hand_side)
- {
+ {
const unsigned int n_q_points
= quadrature.n_quadrature_points;
-
+
cell_residual.resize(n_q_points);
- rhs_values.resize(n_q_points);
- dual_weights.resize(n_q_points);
+ rhs_values.resize(n_q_points);
+ dual_weights.resize(n_q_points);
cell_grad_grads.resize(n_q_points);
}
-
-
+
+
template <int dim>
WeightedResidual<dim>::FaceData::
update_normal_vectors),
fe_subface_values_cell (fe, face_quadrature,
update_gradients)
- {
+ {
const unsigned int n_face_q_points
= face_quadrature.n_quadrature_points;
-
+
jump_residual.resize(n_face_q_points);
- dual_weights.resize(n_face_q_points);
+ dual_weights.resize(n_face_q_points);
cell_grads.resize(n_face_q_points);
neighbor_grads.resize(n_face_q_points);
}
-
+
threads.join_all ();
}
-
+
template <int dim>
void
WeightedResidual<dim>::solve_primal_problem ()
{
DualSolver<dim>::solve_problem ();
}
-
+
template <int dim>
void
{
PrimalSolver<dim>::postprocess (postprocessor);
}
-
-
+
+
template <int dim>
unsigned int
WeightedResidual<dim>::n_dofs () const
0.8, 0.02);
this->triangulation->execute_coarsening_and_refinement ();
}
-
+
template <int dim>
void
{
const PrimalSolver<dim> &primal_solver = *this;
const DualSolver<dim> &dual_solver = *this;
-
+
ConstraintMatrix primal_hanging_node_constraints;
DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
primal_hanging_node_constraints);
dual_solver.solution,
primal_solver.dof_handler,
primal_hanging_node_constraints,
- dual_solution);
+ dual_solution);
DataOut<dim> data_out;
data_out.attach_dof_handler (primal_solver.dof_handler);
"primal_solution");
data_out.add_data_vector (dual_solution,
"dual_solution");
-
+
data_out.build_patches ();
-
+
std::ostringstream filename;
filename << "solution-"
<< this->refinement_cycle
-
+
template <int dim>
void
WeightedResidual<dim>::
dual_solver.dof_handler,
dual_hanging_node_constraints,
primal_solution);
-
+
ConstraintMatrix primal_hanging_node_constraints;
DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
primal_hanging_node_constraints);
primal_solver.dof_handler,
primal_hanging_node_constraints,
dual_weights);
-
-
+
+
FaceIntegrals face_integrals;
for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
cell!=dual_solver.dof_handler.end();
threads.join_all();
}
- unsigned int present_cell=0;
+ unsigned int present_cell=0;
for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
cell!=dual_solver.dof_handler.end();
++cell, ++present_cell)
*dual_solver.quadrature,
*primal_solver.rhs_function);
FaceData face_data (*dual_solver.fe,
- *dual_solver.face_quadrature);
+ *dual_solver.face_quadrature);
active_cell_iterator cell=dual_solver.dof_handler.begin_active();
for (unsigned int t=0;
(t<this_thread) && (cell!=dual_solver.dof_handler.end());
- ++t, ++cell);
+ ++t, ++cell)
+ {
+ }
+
if (cell == dual_solver.dof_handler.end())
return;
-
+
for (unsigned int cell_index=this_thread; true; )
{
integrate_over_cell (cell, cell_index,
dual_weights,
cell_data,
error_indicators);
-
+
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
{
- if (cell->face(face_no)->at_boundary())
+ if (cell->face(face_no)->at_boundary())
{
face_integrals[cell->face(face_no)] = 0;
continue;
};
-
+
if ((cell->neighbor(face_no)->has_children() == false) &&
(cell->neighbor(face_no)->level() == cell->level()) &&
(cell->neighbor(face_no)->index() < cell->index()))
if (cell->at_boundary(face_no) == false)
if (cell->neighbor(face_no)->level() < cell->level())
- continue;
+ continue;
if (cell->face(face_no)->has_children() == false)
primal_solution,
dual_weights,
face_data,
- face_integrals);
+ face_integrals);
else
integrate_over_irregular_face (cell, face_no,
primal_solution,
for (unsigned int t=0;
((t<n_threads) && (cell!=dual_solver.dof_handler.end()));
- ++t, ++cell, ++cell_index);
+ ++t, ++cell, ++cell_index)
+ {
+ }
+
if (cell == dual_solver.dof_handler.end())
break;
};
}
-
+
template <int dim>
void WeightedResidual<dim>::
integrate_over_regular_face (const active_cell_iterator &cell,
const unsigned int
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
const active_cell_iterator neighbor = cell->neighbor(face_no);
- face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
+ face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
face_data.neighbor_grads);
face_data.fe_face_values_cell.get_function_values (dual_weights,
face_data.dual_weights);
-
+
double face_integral = 0;
for (unsigned int p=0; p<n_q_points; ++p)
face_integral += (face_data.jump_residual[p] *
}
-
+
template <int dim>
void WeightedResidual<dim>::
integrate_over_irregular_face (const active_cell_iterator &cell,
n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
const typename DoFHandler<dim>::face_iterator
- face = cell->face(face_no);
+ face = cell->face(face_no);
const typename DoFHandler<dim>::cell_iterator
- neighbor = cell->neighbor(face_no);
+ neighbor = cell->neighbor(face_no);
Assert (neighbor.state() == IteratorState::valid,
ExcInternalError());
Assert (neighbor->has_children(),
const unsigned int
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-
+
for (unsigned int subface_no=0;
subface_no<face->n_children(); ++subface_no)
{
neighbor_neighbor);
face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
face_data.neighbor_grads);
-
+
for (unsigned int p=0; p<n_q_points; ++p)
face_data.jump_residual[p]
= ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
face_data.fe_face_values_neighbor.get_function_values (dual_weights,
face_data.dual_weights);
-
+
double face_integral = 0;
for (unsigned int p=0; p<n_q_points; ++p)
face_integral += (face_data.jump_residual[p] *
double sum = 0;
for (unsigned int subface_no=0;
- subface_no<face->n_children(); ++subface_no)
+ subface_no<face->n_children(); ++subface_no)
{
Assert (face_integrals.find(face->child(subface_no)) !=
face_integrals.end(),
ExcInternalError());
Assert (face_integrals[face->child(subface_no)] != -1e20,
ExcInternalError());
-
+
sum += face_integrals[face->child(subface_no)];
};
face_integrals[face] = sum;
}
-
+
}
typedef std::list<Evaluator*> EvaluatorList;
- struct ProblemDescription
+ struct ProblemDescription
{
unsigned int primal_fe_degree;
unsigned int dual_fe_degree;
*descriptor.dual_functional);
break;
};
-
+
case ProblemDescription::global_refinement:
{
solver
descriptor.data->get_boundary_values());
break;
};
-
+
case ProblemDescription::kelly_indicator:
{
solver
*descriptor.kelly_weight);
break;
};
-
+
default:
AssertThrow (false, ExcInternalError());
};
-
+
for (unsigned int step=0; true; ++step)
{
deallog << "Refinement cycle: " << step
<< std::endl;
-
+
solver->set_refinement_cycle (step);
solver->solve_problem ();
solver->output_solution ();
deallog << " Number of degrees of freedom="
<< solver->n_dofs() << std::endl;
-
+
for (typename EvaluatorList::const_iterator
e = descriptor.evaluator_list.begin();
e != descriptor.evaluator_list.end(); ++e)
solver->postprocess (**e);
};
-
+
if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
solver->refine_grid ();
else
-int main ()
+int main ()
{
deallog.depth_console (0);
try
{
deallog << std::setprecision(2);
logfile << std::setprecision(2);
-
+
deallog.attach(logfile);
deallog.depth_console(0);
- deallog.threshold_double(1.e-10);
+ deallog.threshold_double(1.e-10);
const unsigned int dim = 2;
Framework<dim>::ProblemDescription descriptor;
descriptor.dual_fe_degree = 2;
descriptor.data = new Data::SetUp<Data::Exercise_2_3<dim>,dim> ();
-
+
const Point<dim> evaluation_point (0.75, 0.75);
descriptor.dual_functional
= new DualFunctional::PointValueEvaluation<dim> (evaluation_point);
-
+
Evaluation::PointValueEvaluation<dim>
postprocessor1 (evaluation_point);
Evaluation::GridOutput<dim>
postprocessor2 ("grid");
-
+
descriptor.evaluator_list.push_back (&postprocessor1);
descriptor.evaluator_list.push_back (&postprocessor2);
descriptor.max_degrees_of_freedom = 1000;
-
+
Framework<dim>::run (descriptor);
}
<< std::endl;
return 1;
}
- catch (...)
+ catch (...)
{
deallog << std::endl << std::endl
<< "----------------------------------------------------"