-// deal_II_libraries.g=-ldeal_II_2d.g
-// deal_II_libraries=-ldeal_II_2d
+//---------------------------- wave-test-3.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- wave-test-3.cc ---------------------------
+
#include <grid/tria_boundary_lib.h>
#include <base/parameter_handler.h>
template <int dim> class TimeStep_Postprocess;
-
template <int dim>
class TimeStepBase_Wave : public TimeStepBase_Tria<dim>{
public:
};
-
template <int dim>
class TimeStep_Wave : public virtual TimeStepBase_Wave<dim>
{
virtual string branch_signature () const = 0;
DeclException0 (ExcIO);
DeclException0 (ExcCoarsestGridsDiffer);
-
-
- protected:
+
+
+protected:
struct StatisticData
{
StatisticData ();
};
-
template <int dim>
class TimeStep_Primal : public TimeStep_Wave<dim>
{
};
-
template <int dim>
class TimeStep_Dual : public TimeStep_Wave<dim>
{
};
-
#include <lac/full_matrix.h>
-
template <int dim>
class TimeStep_ErrorEstimation : public virtual TimeStepBase_Wave<dim>
{
};
- struct CellwiseError
+struct CellwiseError
{
CellwiseError (const unsigned int n_errors);
vector<ErrorOnCell> errors;
};
-
template <int dim> class WaveParameters;
-
template <int dim>
class TimeStep : public TimeStep_Primal<dim>, public TimeStep_Dual<dim>, public TimeStep_Postprocess<dim>
{
template <int dim> class TimeStep_Dual;
-
template <int dim>
class DualFunctional {
public:
};
-
-
-
/**
* Compute the dual functional which is approximately associated
* with the end time energy in the high atmosphere above 4000km.
};
-
-
-
-
-
/**
* Let the point value of $u$ at the origin integrated over time
* be the goal.
};
-
-
-
/**
* Dual function corresponding to the #EvaluateSeismicSignal# class.
*/
};
-
/**
* Compute the dual problem associated with the functional
* $J(\Psi) = \int u ds$ with the integral being over some
};
-
/**
* Compute $J(\Psi) = \int_0^0.25 u(x=2,y,t=2.2) dy.
*/
};
-
-
-
/**
* 1d test case, evaluating the region (-.5,.5) at the endtime. Intended for some
* tests on split triangulations with one fine and one coarse region.
};
-
/**
* Compute $J(\Psi) = \int_{-0.6}^{-0.4} u(x,t=2.5) dx.
*/
};
-
/**
* Compute $J(\Psi) = \int_{-0.1}^{0.1} u(x,t=2.4) dx.
*/
};
-
/**
*/
template <int dim>
#include <string>
-
-
-
/**
* This class provides a simple interface to do arbitrary evaluations of
* the numerical solution. Concrete classes implementing evaluations
};
-
-
-
/**
* This class is a common base class to the following two. It provides
* for some infrastructure for evaluations computing the energy in part
};
-
-
-
-
-
/**
* Evaluate the value of $u$ at the origin, i.e. $u(0,0)$.
*
};
-
-
-
-
/**
* Integrate the value of $u_h$ at the top boundary over $x$ and $t$ using a
* highly oscillatory weight.
const double pi = 3.14159265359;
return sin(3*pi*p(0))*sin(5*pi*time/2);
};
-
-
- virtual double evaluate ();
+
+
+virtual double evaluate ();
virtual void print_final_result (ostream &out);
virtual double get_final_result ();
virtual string description () const;
};
-
/**
* Integrate the value of $u_h$ at the top line $x=1.5, y=0..1/16$ at $t=1.6..1.8$.
*/
public:
EvaluateSplitSignal () :
result (0) {};
-
-
- virtual double evaluate ();
+
+
+virtual double evaluate ();
virtual void print_final_result (ostream &out);
virtual double get_final_result ();
virtual string description () const;
double result;
};
-
template <int dim>
class EvaluateOneBranch1d : public EvaluationBase<dim> {
public:
EvaluateOneBranch1d () :
result (0) {};
-
-
- virtual double evaluate ();
+
+
+virtual double evaluate ();
virtual void print_final_result (ostream &out);
virtual double get_final_result ();
virtual string description () const;
};
-
-
template <int dim>
class EvaluateSecondCrossing1d : public EvaluationBase<dim> {
public:
EvaluateSecondCrossing1d () :
result (0) {};
-
-
- virtual double evaluate ();
+
+
+virtual double evaluate ();
virtual void print_final_result (ostream &out);
virtual double get_final_result ();
virtual string description () const;
};
-
template <int dim>
class EvaluateHuyghensWave : public EvaluationBase<dim> {
public:
EvaluateHuyghensWave () :
integrated_value (0),
weighted_value (0) {};
-
-
- virtual double evaluate ();
+
+
+virtual double evaluate ();
virtual void print_final_result (ostream &out);
virtual double get_final_result ();
virtual string description () const;
};
-
-
template <int dim> class DataOutStack;
-
/**
* This class has some data members which are shared between the different
* time steps within one sweep. Unlike the #SweepInfo# class, the members
};
-
#include <base/timer.h>
#include <iostream>
#include <list>
-
/**
* This class provides some data members which collect information on the
* different time steps of one sweep.
};
- struct Timers
+struct Timers
{
Timer grid_generation;
Timer primal_problem;
};
- Data & get_data ();
+Data & get_data ();
Timers & get_timers ();
-
-
- template <int dim>
+
+
+template <int dim>
void write_summary (const list<EvaluationBase<dim>*> & eval_list,
ostream &out) const;
};
-
#include <lac/sparse_matrix.h>
-
-
/**
* Enum denoting the different possibilities to precondition a solver.
*/
};
-
/**
* Wrapper for the #SparseMatrix<double># class which handles the preconditioning.
*/
};
-
/**
* This is a helper class which has a collection of static elements and returns
* the right finite element as a pointer when the name of the element is given.
};
-
#include <base/forward_declarations.h>
#include <grid/forward_declarations.h>
#include <list>
template <int dim> class EvaluationBase;
-
-
/**
* This is a class holding all the input parameters to the program. It is more
* or less a loose collection of data and the only purpose of this class is
* Undefined copy operator.
*/
WaveParameters & operator = (const WaveParameters &);
-
-
- /**
+
+
+/**
* List of names for the initial values.
* Make this a member of the templated
* class since the supported initial
*/
static const string dual_functional_names;
-
- /**
+
+/**
* Set the initial function pointers
* depending on the given names.
*/
};
-
#include <numerics/time_dependent.h>
template <int dim> class WaveParameters;
class SweepInfo;
-
/**
* Top-level class of the timestepping mechanism. This class manages
* the execution and solution of primal and dual problem, of computing
};
-
-
-
/**
* Top-level class providing the set up of a simulation. The
* class provides an interface suitable to the #MultipleParameterLoop#
};
-
-
-
/* $Id$ */
#include <cmath>
-
/*------------------------ DualFunctional --------------------------------*/
template <int dim>
{};
-
template <int dim>
void DualFunctional<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
template <int dim>
void DualFunctional<dim>::compute_endtime_vectors (Vector<double> &final_u_bar,
Vector<double> &final_v_bar) {
};
-
template <int dim>
bool DualFunctional<dim>::use_primal_solutions () const {
return use_primal_problem;
};
-
template <int dim>
bool DualFunctional<dim>::use_primal_solutions_at_endtime () const {
return use_primal_problem_at_endtime;
};
-
template <int dim>
void DualFunctional<dim>::reset (const TimeStep_Primal<dim> &primal_problem) {
Assert (use_primal_problem ||
};
-
template <int dim>
void DualFunctional<dim>::reset (const TimeStep_Dual<dim> &dual_problem) {
tria = dual_problem.tria;
};
-
-
-
-
-
-
-
/* ----------------------- EndEnergy ------------------------------*/
DualFunctional<dim> (use_primal_problem, true) {};
-
-
template <int dim>
void EndEnergy<dim>::compute_vectors (const PartOfDomain pod,
Vector<double> &final_u_bar,
continue;
break;
};
-
-
- fe_values.reinit (cell);
+
+
+fe_values.reinit (cell);
fe_values_primal.reinit (primal_cell);
fe_values_primal.get_function_values (*v, local_v);
fe_values_primal.get_function_grads (*u, local_u_grad);
};
-
-
-
/*------------------------ IntegrateValueAtOrigin --------------------------------*/
};
-
-
-
-
/*------------------------ SeismicSignal --------------------------------*/
};
-
template <int dim>
void SeismicSignal<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
-
/*------------------------ EarthSurface --------------------------------*/
};
-
template <int dim>
void EarthSurface<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
-
/*------------------------ SplitSignal --------------------------------*/
};
-
-
template <int dim>
void SplitSignal<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
-
/* ------------------------------ Split line 1d case ----------------------------- */
template <int dim>
#endif
-
/*------------------------ OneBranch1d --------------------------------*/
-
template <int dim>
void OneBranch1d<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
/*------------------------ SecondCrossing --------------------------------*/
-
template <int dim>
void SecondCrossing<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
/*------------------------ HuyghensWave --------------------------------*/
-
template <int dim>
void HuyghensWave<dim>::compute_functionals (Vector<double> &j1,
Vector<double> &j2) {
};
-
// explicit specializations
template class DualFunctional<2>;
#include <fstream>
-
/*--------------------------- EvaluationBase --------------------------*/
template <int dim>
};
-
template <int dim>
void EvaluationBase<dim>::reset () {};
-
template <int dim>
void EvaluationBase<dim>::print_final_result (ostream &) {};
};
-
-
-
-
/*--------------------------- EvaluateEnergyContent ----------------------*/
template <int dim>
};
-
template <int dim>
double EvaluateEnergyContent<dim>::compute_energy (const PartOfDomain pod) const {
const double y_offset = 300000000;
continue;
break;
};
-
-
- fe_values.reinit (cell);
+
+
+fe_values.reinit (cell);
const FullMatrix<double> &values = fe_values.get_shape_values();
const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
const vector<double> &weights = fe_values.get_JxW_values ();
};
-
-
/* ---------------------------- EvaluateIntegratedValueAtOrigin ------------------- */
};
-
template <int dim>
double EvaluateIntegratedValueAtOrigin<dim>::get_final_result () {
return integrated_value;
};
-
template <int dim>
string EvaluateIntegratedValueAtOrigin<dim>::description () const {
return "integrated value at origin";
};
-
template <int dim>
void EvaluateIntegratedValueAtOrigin<dim>::reset () {
integrated_value = 0;
};
-
template <int dim>
double EvaluateIntegratedValueAtOrigin<dim>::evaluate () {
DoFHandler<dim>::active_cell_iterator cell = dof->begin_active(),
};
-
-
-
-
/*------------------------- EvaluateSeismicSignal --------------------------*/
};
-
template <int dim>
double EvaluateSeismicSignal<dim>::get_final_result () {
return result;
};
-
template <int dim>
string EvaluateSeismicSignal<dim>::description () const {
return "Integrated seismic signal at top";
};
-
template <int dim>
void EvaluateSeismicSignal<dim>::reset () {
result = 0;
};
-
template <>
double EvaluateSeismicSignal<1>::evaluate ()
{
};
-
template <int dim>
double EvaluateSeismicSignal<dim>::evaluate () {
const unsigned int n_q_points = quadrature_face->n_quadrature_points;
};
-
-
/*------------------------- EvaluateSplitSignal --------------------------*/
};
-
template <int dim>
double EvaluateSplitSignal<dim>::get_final_result () {
return result;
};
-
template <int dim>
string EvaluateSplitSignal<dim>::description () const {
return "Integrated split signal (exact: (2+pi)/(16-pi)=0.010229)";
};
-
template <int dim>
void EvaluateSplitSignal<dim>::reset () {
result = 0;
};
-
template <>
double EvaluateSplitSignal<1>::evaluate ()
{
};
-
template <int dim>
double EvaluateSplitSignal<dim>::evaluate () {
if ((time<=1.6) || (time>1.8))
};
-
-
/*------------------------- EvaluateOneBranch1d --------------------------*/
};
-
template <int dim>
double EvaluateOneBranch1d<dim>::get_final_result () {
return result;
};
-
template <int dim>
string EvaluateOneBranch1d<dim>::description () const {
return "One branch integrated (exact: 0.055735)";
};
-
template <int dim>
void EvaluateOneBranch1d<dim>::reset () {
result = 0;
};
-
template <int dim>
double EvaluateOneBranch1d<dim>::evaluate ()
{
#endif
-
-
/*------------------------- EvaluateSecondCrossing1d --------------------------*/
};
-
template <int dim>
double EvaluateSecondCrossing1d<dim>::get_final_result () {
return result;
};
-
template <int dim>
string EvaluateSecondCrossing1d<dim>::description () const {
return "Second crossing (exact: 0.011147)";
};
-
template <int dim>
void EvaluateSecondCrossing1d<dim>::reset () {
result = 0;
};
-
template <int dim>
double EvaluateSecondCrossing1d<dim>::evaluate ()
{
#endif
-
/*------------------------- EvaluateHuyghensWave --------------------------*/
};
-
template <int dim>
double EvaluateHuyghensWave<dim>::get_final_result () {
return weighted_value / integrated_value;
};
-
template <int dim>
string EvaluateHuyghensWave<dim>::description () const {
return "Huyghens wave";
};
-
template <int dim>
void EvaluateHuyghensWave<dim>::reset () {
integrated_value = weighted_value = 0;
};
-
template <int dim>
double EvaluateHuyghensWave<dim>::evaluate ()
{
};
-
// explicit instantiations
template class EvaluationBase<2>;
template class EvaluateEnergyContent<2>;
{};
-
template <int dim>
void TimestepManager<dim>::run_sweep (const unsigned int sweep_no)
{
(parameters.refinement_strategy == WaveParameters<dim>::dual_estimator))
sweep_data.data_out_stack->declare_data_vector ("est_error", DataOutStack<dim>::cell_vector);
};
-
-
- deallog << "Sweep " << setw(2) << sweep_no << ':' << endl
+
+
+deallog << "Sweep " << setw(2) << sweep_no << ':' << endl
<< "---------" << endl;
for (typename list<EvaluationBase<dim>*>::const_iterator i = parameters.eval_list.begin();
};
-
-
template <int dim>
void TimestepManager<dim>::refine_grids ()
{
deallog << " Collecting refinement data: " << endl;
-
- const unsigned int n_timesteps = timesteps.size();
+
+const unsigned int n_timesteps = timesteps.size();
// first collect all the error indicators
vector<Vector<float> > indicators (n_timesteps);
->get_timestep_postprocess().get_tria_refinement_criteria (indicators[i]);
- // count the number of cells for some
+// count the number of cells for some
// statistics and other things
unsigned int total_number_of_cells = 0;
for (unsigned int i=0; i<timesteps.size(); ++i)
total_number_of_cells += indicators[i].size();
-
- if (parameters.produce_error_statistics)
+if (parameters.produce_error_statistics)
{
deallog << " Generating error statistics ";
};
- if (parameters.compare_indicators_globally)
+if (parameters.compare_indicators_globally)
{
// collect all indicators in one
static_cast<TimeStepBase_Tria<dim>*>(timesteps.back())->sleep(1);
- deallog << " Got " << total_number_of_cells << " presently, expecting "
+deallog << " Got " << total_number_of_cells << " presently, expecting "
<< total_expected_cells << " for next sweep." << endl;
};
};
-
-
template <int dim>
void TimestepManager<dim>::write_statistics (const SweepInfo &sweep_info) const
{
deallog << endl;
};
-
- // write summary
+
+// write summary
if (true)
{
deallog << " Writing summary.";
};
};
-
template <int dim>
void TimestepManager<dim>::write_stacked_data (DataOutStack<dim> &data_out_stack) const
};
-
-
//explicit instantiation
template class TimestepManager<2>;
"|Huyghens wave");
-
DeclException1 (ExcUnknownName,
string,
<< "Unknown description string " << arg1);
-
template <int dim>
class InitialValues {
public:
};
- class SmallBump : public Function<dim> {
+class SmallBump : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
- class ShiftedBump : public Function<dim> {
+class ShiftedBump : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
- class Earthquake : public Function<dim> {
+class Earthquake : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
-
template <int dim>
class Coefficients {
public:
};
- class Gradient : public Function<dim> {
+class Gradient : public Function<dim> {
public:
inline virtual double value (const Point<dim> &p,
const unsigned int) const {
};
-
- class PreliminaryEarthModel : public Function<dim> {
+class PreliminaryEarthModel : public Function<dim> {
public:
inline virtual double value (const Point<dim> &p,
const unsigned int) const {
};
};
-
-
- class Distorted : public Function<dim> {
+
+class Distorted : public Function<dim> {
public:
inline virtual double value (const Point<dim> &p,
const unsigned int) const {
};
-
-
template <int dim>
class BoundaryValues {
public:
};
- class FastWaveFromLeft_u : public Function<dim> {
+class FastWaveFromLeft_u : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
};
-
- class WaveFromLeftCenter_u : public Function<dim> {
+
+class WaveFromLeftCenter_u : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
- class WaveFromLeftBottom_u : public Function<dim> {
+class WaveFromLeftBottom_u : public Function<dim> {
public:
virtual double value (const Point<dim> &p,
const unsigned int) const {
};
-
template <int dim>
class Boundaries
{
middle *= sqrt(line->vertex(0).square()) / sqrt(middle.square());
return middle;
};
-
-
- virtual Point<dim>
+
+
+virtual Point<dim>
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const {
Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
middle *= sqrt(quad->vertex(0).square()) / sqrt(middle.square());
};
-
-
template <int dim>
WaveParameters<dim>::WaveParameters () :
boundary_values_u (0),
{};
-
template <int dim>
WaveParameters<dim>::~WaveParameters ()
{
};
-
template <int dim>
void WaveParameters<dim>::delete_parameters ()
{
};
-
-
template <int dim>
void WaveParameters<dim>::set_initial_functions (const string &u_name,
const string &v_name) {
};
-
-
-
template <int dim>
void WaveParameters<dim>::set_coefficient_functions (const string &name) {
Assert (density==0, ExcInternalError());
};
-
template <int dim>
void WaveParameters<dim>::set_boundary_functions (const string &name) {
Assert (boundary_values_u==0, ExcInternalError());
};
-
template <int dim>
void WaveParameters<dim>::make_eval_list (const string &names) {
Assert (eval_list.size()==0, ExcInternalError());
};
-
-
template <int dim>
void WaveParameters<dim>::set_dual_functional (const string &name) {
Assert (dual_functional==0, ExcInternalError());
};
-
-
#if 2 == 1
template <>
#endif
-
#if 2 == 2
template <>
break;
};
-
-
- case uniform_channel:
+
+
+case uniform_channel:
{
coarse_grid->refine_global (initial_refinement-1);
break;
};
- default:
+default:
Assert (false, ExcInternalError());
};
};
break;
};
-
- case square:
+
+case square:
case seismic_square:
{
GridGenerator::hyper_cube (*coarse_grid, -1, 1);
#endif
-
-
template <int dim>
void WaveParameters<dim>::declare_parameters (ParameterHandler &prm)
{
prm.leave_subsection ();
- prm.enter_subsection ("Goal");
+prm.enter_subsection ("Goal");
prm.declare_entry ("Goal", "none",
Patterns::Selection (dual_functional_names));
prm.declare_entry ("Evaluate", "");
prm.leave_subsection ();
- prm.declare_entry ("Refinement criterion", "energy estimator",
+prm.declare_entry ("Refinement criterion", "energy estimator",
Patterns::Selection ("energy estimator|dual estimator"));
prm.declare_entry ("Sweeps", "3", Patterns::Integer());
};
-
template <int dim>
void WaveParameters<dim>::parse_parameters (ParameterHandler &prm) {
// declare some maps for convenience,
write_strategy_list["never"] = never;
write_strategy_list["all sweeps"] = all_sweeps;
write_strategy_list["last sweep only"] = last_sweep_only;
-
-
- prm.enter_subsection ("Grid");
+
+
+prm.enter_subsection ("Grid");
initial_refinement = prm.get_integer ("Initial refinement");
// don't make the grid here already, since
// it may depend on the chosen boundary
prm.leave_subsection ();
- prm.enter_subsection ("Goal");
+prm.enter_subsection ("Goal");
set_dual_functional (prm.get("Goal"));
make_eval_list (prm.get("Evaluate"));
prm.leave_subsection ();
-
-
- if (prm.get("Refinement criterion")=="energy estimator")
+
+if (prm.get("Refinement criterion")=="energy estimator")
refinement_strategy = energy_estimator;
else
refinement_strategy = dual_estimator;
};
-
-
// explicit instantiations
template class WaveParameters<2>;
/* $Id$ */
};
-
template <int dim>
SweepData<dim>::~SweepData ()
{
};
-
-
// explicit instantiations
template class SweepData<2>;
/* $Id$ */
};
-
SweepInfo::Timers &
SweepInfo::get_timers ()
{
};
-
template <int dim>
void
SweepInfo::write_summary (const list<EvaluationBase<dim>*> &eval_list,
};
-
-
-
-
-
SweepInfo::Data::Data () :
accumulated_error (0),
cells (0),
{};
-
-
-
-
-
// explicit instantiations
template
void SweepInfo::write_summary (const list<EvaluationBase<2>*> &eval_list,
#include <iomanip>
-
-
static const pair<unsigned int, double> relaxations[3]
= { make_pair(100,5), make_pair(300,3), make_pair(500,2) };
&relaxations[3]));
-
template <int dim>
TimeStepBase_Wave<dim>::TimeStepBase_Wave ():
TimeStepBase_Tria<dim> (),
{};
-
template <int dim>
TimeStepBase_Wave<dim>::TimeStepBase_Wave (const double time,
TimeStepBase_Tria<dim>::Flags flags,
{};
-
template <int dim>
const TimeStep_Primal<dim> &
TimeStepBase_Wave<dim>::get_timestep_primal () const
};
-
template <int dim>
const TimeStep_Dual<dim> &
TimeStepBase_Wave<dim>::get_timestep_dual () const
};
-
template <int dim>
const TimeStep_Postprocess<dim> &
TimeStepBase_Wave<dim>::get_timestep_postprocess () const
};
-
template <int dim>
string TimeStepBase_Wave<dim>::tmp_filename_base (const string &branch_signature) const
{
};
-
template <int dim>
void TimeStepBase_Wave<dim>::attach_sweep_info (SweepInfo &si)
{
};
-
template <int dim>
void TimeStepBase_Wave<dim>::attach_sweep_data (SweepData<dim> &sd)
{
};
-
-
-
-
/* --------------------------------------------------------------*/
{};
-
template <int dim>
TimeStep_Wave<dim>::~TimeStep_Wave ()
{
Assert (v.size() ==0, ExcInternalError());
};
-
template <int dim>
void TimeStep_Wave<dim>::wake_up (const unsigned int wakeup_level)
if (parameters.renumber_dofs)
DoFRenumbering::Cuthill_McKee (*dof_handler);
-
- constraints.clear ();
+
+constraints.clear ();
DoFTools::make_hanging_node_constraints (*dof_handler, constraints);
constraints.close ();
};
-
template <int dim>
void TimeStep_Wave<dim>::sleep (const unsigned int sleep_level)
{
};
-
template <int dim>
void TimeStep_Wave<dim>::end_sweep ()
{
};
-
template <int dim>
unsigned int TimeStep_Wave<dim>::solve (const UserMatrix &matrix,
Vector<double> &solution,
};
-
template <int dim>
void TimeStep_Wave<dim>::create_matrices ()
{
if (stiffness_constant)
fill_n (stiffness_values.begin(), n_q_points,
parameters.stiffness->value(Point<dim>()));
-
-
- FEValues<dim> fe_values (fe, quadrature,
+
+
+FEValues<dim> fe_values (fe, quadrature,
UpdateFlags(update_values |
update_gradients |
update_JxW_values |
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_laplace_matrix (dofs_per_cell, dofs_per_cell);
-
- for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler->begin_active();
+
+for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler->begin_active();
cell != dof_handler->end(); ++cell)
{
fe_values.reinit (cell);
};
-
template <int dim>
void TimeStep_Wave<dim>::transfer_old_solutions (Vector<double> &old_u,
Vector<double> &old_v) const
};
-
template <int dim>
void
TimeStep_Wave<dim>::transfer_old_solutions (const typename DoFHandler<dim>::cell_iterator &old_cell,
};
-
template <int dim>
pair<double,double>
TimeStep_Wave<dim>::compute_energy () {
};
-
template <int dim>
TimeStep_Wave<dim>::StatisticData::
StatisticData () :
{};
-
template <int dim>
TimeStep_Wave<dim>::StatisticData::
StatisticData (const unsigned int n_active_cells,
{};
-
template <int dim>
void
TimeStep_Wave<dim>::StatisticData::write_descriptions (ostream &out)
};
-
template <int dim>
void TimeStep_Wave<dim>::StatisticData::write (ostream &out) const
{
};
-
-
-
-
-
// explicit instantiations
template class TimeStepBase_Wave<2>;
template class TimeStep_Wave<2>;
#include <iomanip>
-
template <int dim>
TimeStep_Dual<dim>::TimeStep_Dual (const string &dual_fe)
:
{};
-
template <int dim>
void TimeStep_Dual<dim>::do_initial_step () {
deallog << " Dual problem: time="
};
-
template <int dim>
void TimeStep_Dual<dim>::do_timestep ()
{
// zero already.
- // in 1d, do not set boundary conditions
+// in 1d, do not set boundary conditions
// at all
//
// note: in boundary_value_map, all entries
};
-
-
template <int dim>
void TimeStep_Dual<dim>::solve_dual_problem ()
{
};
-
template <int dim>
string TimeStep_Dual<dim>::branch_signature () const
{
};
-
template <int dim>
void TimeStep_Dual<dim>::wake_up (const unsigned int wakeup_level)
{
};
-
template <int dim>
void TimeStep_Dual<dim>::assemble_vectors (Vector<double> &right_hand_side1,
Vector<double> &right_hand_side2) {
};
-
template <int dim>
void TimeStep_Dual<dim>::build_rhs (Vector<double> &right_hand_side1,
Vector<double> &right_hand_side2) {
update_JxW_values |
update_q_points));
-
- cell_iterator old_cell = previous_time_level.dof_handler->begin(),
+
+cell_iterator old_cell = previous_time_level.dof_handler->begin(),
new_cell = dof_handler->begin(),
end_cell = (tria->n_levels() == 1 ?
static_cast<cell_iterator>(dof_handler->end()) :
};
-
template <int dim>
void
TimeStep_Dual<dim>::build_rhs (const DoFHandler<dim>::cell_iterator &old_cell,
};
- // select the TimeStep_Wave part in the
+// select the TimeStep_Wave part in the
// TimeStep_Dual branch
const TimeStep_Dual<dim> &previous_time_level
= static_cast<const TimeStepBase_Wave<dim>*>(next_timestep)->get_timestep_dual();
};
-
template <int dim>
unsigned int
TimeStep_Dual<dim>::collect_from_children (const DoFHandler<dim>::cell_iterator &old_cell,
};
-
template <int dim>
unsigned int
TimeStep_Dual<dim>::distribute_to_children (const DoFHandler<dim>::cell_iterator &new_cell,
// the new grid
vector<unsigned int> new_dof_indices (dofs_per_cell, DoFHandler<dim>::invalid_dof_index);
-
- // loop over the child cells
+
+// loop over the child cells
for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
{
const DoFHandler<dim>::cell_iterator new_child = new_cell->child(c);
};
-
-
// explicit instantiations
template class TimeStep_Dual<2>;
/* $Id$ */
#include <numeric>
-
template <int dim>
TimeStep_ErrorEstimation<dim>::TimeStep_ErrorEstimation ()
{};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::estimate_error ()
{
};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::wake_up (const unsigned int wakeup_level)
{
};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::sleep (const unsigned int sleep_level)
{
};
-
template <int dim>
void
TimeStep_ErrorEstimation<dim>::get_tria_refinement_criteria (Vector<float> &indicators) const
};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::estimate_error_energy (const unsigned int which_variables) {
Assert (which_variables<=1, ExcInternalError());
};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::estimate_error_dual () {
CellwiseError cellwise_error (tria->n_active_cells());
(previous_timestep)->get_timestep_dual();
- // first clear the user pointers of
+// first clear the user pointers of
// the cells we need
if (true)
{
// compute the sum of the errors
// on the cells
ErrorOnCell total_estimated_error;
-
- // now fill the data we collected to the
+
+// now fill the data we collected to the
// error_per_cell array
Vector<float>::iterator i = estimated_error_per_cell.begin();
DoFHandler<dim>::active_cell_iterator
};
-
template <int dim>
void
TimeStep_ErrorEstimation<dim>::estimate_error_dual (const DoFHandler<dim>::cell_iterator &primal_cell,
return;
};
-
- const TimeStep_Primal<dim> &primal_problem = get_timestep_primal(),
+const TimeStep_Primal<dim> &primal_problem = get_timestep_primal(),
&primal_problem_old = static_cast<const TimeStepBase_Wave<dim>*>
(previous_timestep)->get_timestep_primal();
const TimeStep_Dual<dim> &dual_problem = get_timestep_dual(),
dofs_per_cell_dual = dual_fe.dofs_per_cell;
- // none of the two cells has children
+// none of the two cells has children
if (!primal_cell->has_children() && !primal_cell_old->has_children())
{
// vector holding the solutions on
dual_cell->get_dof_values (dual_problem.v, local_v_bar);
- // fill local old solution vectors.
+// fill local old solution vectors.
// no problems here, since the two
// cells are both unrefined
primal_cell_old->get_dof_values (primal_problem_old.u, primal_tmp);
};
-
- // only new cell has children. handle this
+// only new cell has children. handle this
// case by prolonging the solutions on the
// old cell to its children and recursing
// thereon
dual_cell_old->get_dof_values (dual_problem_old.v, local_v_bar_old);
- compute_error_on_new_children (primal_cell, dual_cell,
+compute_error_on_new_children (primal_cell, dual_cell,
local_u_old,
local_v_old,
local_u_bar_old,
};
- // last possibility: new cell is not
+// last possibility: new cell is not
// refined, but old one is. in this case:
// collect error on this cell from the
// smaller ones on the old grid
};
- Assert (false, ExcInternalError());
+Assert (false, ExcInternalError());
};
-
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::
compute_error_on_new_children (const DoFHandler<dim>::cell_iterator &primal_cell,
const FiniteElement<dim> &dual_fe = get_timestep_dual().fe;
const unsigned int dofs_per_cell_dual = dual_fe.dofs_per_cell;
-
-
- for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+
+
+for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
{
// we have the solutions on the
// old (large) cell, we restrict it to
};
-
template <int dim>
typename TimeStep_ErrorEstimation<dim>::ErrorOnCell
TimeStep_ErrorEstimation<dim>::collect_error_from_children (const DoFHandler<dim>::cell_iterator &primal_cell_old,
local_difference_u_bar_old -= local_Ih_u_bar_old;
local_difference_v_bar_old = local_v_bar_old;
local_difference_v_bar_old -= local_Ih_v_bar_old;
-
-
- error_sum += error_formula (old_dual_child,
+
+
+error_sum += error_formula (old_dual_child,
child_u, child_v,
child_u_bar, child_v_bar,
local_u_old, local_v_old,
};
-
template <int dim>
typename TimeStep_ErrorEstimation<dim>::ErrorOnCell
TimeStep_ErrorEstimation<dim>::error_formula (const DoFHandler<dim>::active_cell_iterator &cell,
};
-
template <int dim>
typename TimeStep_ErrorEstimation<dim>::ErrorOnCell
TimeStep_ErrorEstimation<dim>::error_formula (const DoFHandler<dim>::active_cell_iterator &cell,
Vector<double> tmp2(dofs_per_cell);
- vector<double> stiffness(fe_values.n_quadrature_points);
+vector<double> stiffness(fe_values.n_quadrature_points);
parameters.stiffness->value_list (fe_values.get_quadrature_points(),
stiffness);
vector<Tensor<1,dim> > grad_stiffness(fe_values.n_quadrature_points);
weights[point] *
stiffness[point];
};
-
-
-
- // ////////////////////////////////////
+
+
+// ////////////////////////////////////
// Compute the different contributions
// separately. Note that the parts
// 1 and 2a+2b together should give
// a small quantity, since they form
// the first domain residual, which is
// small for elements of odd order.
-
-
- // ////////////////////////////////////
+
+
+// ////////////////////////////////////
// PART 1
//
// let #tmp_dual2# hold the contribution
tmp1 -= local_u_old;
error_on_cell.part[0] = mass_matrix.matrix_scalar_product (tmp1, tmp2);
-
-
- // same thing for the second part
+
+
+// same thing for the second part
// with v instead of u
tmp2 = local_difference_v_bar;
tmp2 += local_difference_v_bar_old;
tmp1 -= local_v_old;
error_on_cell.part[1] = mass_matrix.matrix_scalar_product (tmp1, tmp2);
-
-
- // ////////////////////////////////
+
+// ////////////////////////////////
// PART 2a
//
// let tmp2=(1-I)(u_bar^n+u_bar^(n-1))
error_on_cell.part[2] = -(get_backward_timestep() / 4 *
mass_matrix.matrix_scalar_product (tmp1, tmp2));
-
- // ////////////////////////////////
+
+// ////////////////////////////////
// PART 2b
//
// let tmp1 = v^n-v^(n-1)
mass_matrix.matrix_scalar_product (tmp1, tmp2));
- // ////////////////////////////////
+// ////////////////////////////////
// PART 3a
//
// let tmp2=(1-I)(v_bar^n+v_bar^(n-1))
laplace_matrix.matrix_scalar_product (tmp1, tmp2));
- // ////////////////////////////////
+// ////////////////////////////////
// PART 3b
//
// let tmp1 = u^n-u^(n-1)
laplace_matrix.matrix_scalar_product (tmp1, tmp2));
-
// // ///////////////////////////
// // PART 0:
// // tmp1 = u^n-u^(n-1)
// error_on_cell.part[3] = mass_matrix.matrix_scalar_product (tmp1, tmp2);
// error_on_cell.part[3] *= get_backward_timestep() / 12;
-
+
// // /////////////////////////////
// // PART 4 and 6:
// error_on_cell.part[7] = grad_grad_matrix.matrix_scalar_product (tmp1, tmp2);
// error_on_cell.part[7] *= get_backward_timestep() / 12;
-
- return error_on_cell;
-};
-
-
-
+return error_on_cell;
+};
template <int dim>
embedding_matrix(i,j) = primal_fe.shape_value (j, unit_support_points[i]);
- // construct the difference between the
+// construct the difference between the
// identity and the interpolation operator
// to the primal ansatz space. The
// interpolation operator is to act from
};
-
-
-
template <int dim>
TimeStep_ErrorEstimation<dim>::StatisticData::StatisticData () :
estimated_error (0)
{};
-
template <int dim>
TimeStep_ErrorEstimation<dim>::StatisticData::StatisticData (const double estimated_error) :
estimated_error (estimated_error)
{};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::StatisticData::write_descriptions (ostream &out)
{
};
-
template <int dim>
void TimeStep_ErrorEstimation<dim>::StatisticData::write (ostream &out) const
{
};
-
template <int dim>
TimeStep_ErrorEstimation<dim>::ErrorOnCell::ErrorOnCell () {
for (unsigned int i=0; i<sizeof(part)/sizeof(part[0]); ++i)
};
-
template <int dim>
typename TimeStep_ErrorEstimation<dim>::ErrorOnCell
TimeStep_ErrorEstimation<dim>::ErrorOnCell::operator += (const ErrorOnCell &eoc) {
};
-
template <int dim>
double TimeStep_ErrorEstimation<dim>::ErrorOnCell::sum () const {
double x=0;
};
-
template <int dim>
TimeStep_ErrorEstimation<dim>::CellwiseError::CellwiseError (const unsigned int n_errors) :
errors (n_errors),
{};
-
// explicit instantiations
template class TimeStep_ErrorEstimation<2>;
/* $Id$ */
#include <dofs/dof_handler.h>
-
-
template <int dim>
TimeStep<dim>::TimeStep (const double time,
const WaveParameters<dim> ¶meters):
{};
-
-
template <int dim>
void TimeStep<dim>::wake_up (const unsigned int wakeup_level)
{
};
-
template <int dim>
void TimeStep<dim>::sleep (const unsigned int sleep_level)
{
};
-
template <int dim>
void TimeStep<dim>::end_sweep ()
{
};
-
template <int dim>
void TimeStep<dim>::write_statistics_descriptions (ostream &out,
const WaveParameters<dim> ¶meters)
};
-
template <int dim>
void TimeStep<dim>::write_statistics (ostream &out) const
{
};
-
// explicit instantiations
template class TimeStep<2>;
/* $Id$ */
#include <iomanip>
-
template <int dim>
void TimeStep_Postprocess<dim>::postprocess_timestep ()
{
sweep_data->data_out_stack->build_patches ();
sweep_data->data_out_stack->finish_parameter_value ();
};
-
-
- deallog << endl;
+
+
+deallog << endl;
sweep_info->get_timers().postprocessing.stop();
};
-
template <int dim>
void TimeStep_Postprocess<dim>::wake_up (const unsigned int wakeup_level)
{
};
-
template <int dim>
void TimeStep_Postprocess<dim>::sleep (const unsigned int sleep_level)
{
};
-
template <int dim>
string TimeStep_Postprocess<dim>::branch_signature () const
{
};
-
template <int dim>
void TimeStep_Postprocess<dim>::end_sweep ()
{
};
-
-
-
-
template <int dim>
void TimeStep_Postprocess<dim>::interpolate_dual_solution (Vector<double> &interpolated_u_bar,
Vector<double> &interpolated_v_bar) const {
};
-
-
template <int dim>
void TimeStep_Postprocess<dim>::StatisticData::
write_descriptions (ostream &out,
};
-
template <int dim>
void TimeStep_Postprocess<dim>::StatisticData::write (ostream &out) const
{
};
-
-
// explicit instantiations
template class TimeStep_Postprocess<2>;
/* $Id$ */
#include <iomanip>
-
template <int dim>
TimeStep_Primal<dim>::TimeStep_Primal (const string &primal_fe)
:
TimeStep_Wave<dim> (primal_fe)
{};
-
-
template <int dim>
void TimeStep_Primal<dim>::do_initial_step ()
<< dof_handler->n_dofs() << " dofs";
- // add up sweep-accumulated data. count
+// add up sweep-accumulated data. count
// u and v as separate dofs
sweep_info->get_data().cells += tria->n_active_cells();
sweep_info->get_data().primal_dofs += dof_handler->n_dofs() * 2;
};
-
-
template <int dim>
void TimeStep_Primal<dim>::do_timestep ()
{
// u and v as separate dofs
sweep_info->get_data().cells += tria->n_active_cells();
sweep_info->get_data().primal_dofs += dof_handler->n_dofs() * 2;
-
-
- const double time_step = get_backward_timestep ();
+
+
+const double time_step = get_backward_timestep ();
// Vectors holding the right hand sides of
// the two equations.
transfer_old_solutions (old_u, old_v);
};
-
- assemble_vectors (right_hand_side1, right_hand_side2);
+
+assemble_vectors (right_hand_side1, right_hand_side2);
UserMatrix system_matrix (system_sparsity, parameters.preconditioning);
system_matrix.copy_from (mass_matrix);
constraints.condense (right_hand_side2);
- // in 1d, do not set boundary conditions
+// in 1d, do not set boundary conditions
// at all
if (dim != 1)
{
system_matrix, v,
right_hand_side2);
};
-
- if (parameters.extrapolate_old_solutions)
+
+if (parameters.extrapolate_old_solutions)
// solve with a hopefully good guess
// as start vector
{
};
-
-
template <int dim>
void TimeStep_Primal<dim>::solve_primal_problem ()
{
};
-
template <int dim>
string TimeStep_Primal<dim>::branch_signature () const
{
};
-
template <int dim>
void TimeStep_Primal<dim>::wake_up (const unsigned int wakeup_level)
{
};
-
template <int dim>
void TimeStep_Primal<dim>::assemble_vectors (Vector<double> &right_hand_side1,
Vector<double> &right_hand_side2) {
};
-
-
-
-
-
-
template <int dim>
void TimeStep_Primal<dim>::build_rhs (Vector<double> &right_hand_side1,
Vector<double> &right_hand_side2) {
update_JxW_values |
update_q_points));
-
- cell_iterator old_cell = previous_time_level.dof_handler->begin(),
+
+cell_iterator old_cell = previous_time_level.dof_handler->begin(),
new_cell = dof_handler->begin(),
end_cell = (tria->n_levels() == 1 ?
static_cast<cell_iterator>(dof_handler->end()) :
};
-
template <int dim>
void
TimeStep_Primal<dim>::build_rhs (const DoFHandler<dim>::cell_iterator &old_cell,
};
- // select the TimeStep_Wave part in the
+// select the TimeStep_Wave part in the
// TimeStep_Primal branch
const TimeStep_Primal<dim> &previous_time_level
= static_cast<const TimeStepBase_Wave<dim>*>(previous_timestep)->get_timestep_primal();
cell_matrix.vmult (local_A_u, old_dof_values_u);
- rhs1 = local_M_u;
+rhs1 = local_M_u;
rhs1.add (time_step, local_M_v);
rhs1.add ((-time_step*time_step*
parameters.theta*
};
-
-
template <int dim>
unsigned int
TimeStep_Primal<dim>::collect_from_children (const DoFHandler<dim>::cell_iterator &old_cell,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const double time_step = get_backward_timestep();
-
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+
+FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
// these will hold the values of the
// solution on the old grid, i.e. on
};
-
-
template <int dim>
unsigned int
TimeStep_Primal<dim>::distribute_to_children (const DoFHandler<dim>::cell_iterator &new_cell,
// the new grid
vector<unsigned int> new_dof_indices (dofs_per_cell, DoFHandler<dim>::invalid_dof_index);
-
- // loop over the child cells
+
+// loop over the child cells
for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
{
const DoFHandler<dim>::cell_iterator new_child = new_cell->child(c);
};
-
-
-
-
-
// explicit instantiations
template class TimeStep_Primal<2>;
/* $Id$ */
#include <lac/vector.h>
-
void UserMatrix::precondition (Vector<double> &dst,
const Vector<double> &src) const {
switch (preconditioning)
/* $Id$ */
-
#include <fe/fe_lib.lagrange.h>
#include <base/quadrature_lib.h>
-
-
-
-
-
// static objects
const FEQ1<2> FEHelper<2>::fe_linear;
};
-
template <int dim>
const Quadrature<dim> &FEHelper<dim>::get_quadrature (const string &name) {
if (name=="linear")
};
-
template <>
const Quadrature<0> &FEHelper<1>::get_quadrature_face (const string &) {
static const Quadrature<0> dummy_quadrature(1);
};
-
template <int dim>
const Quadrature<dim-1> &FEHelper<dim>::get_quadrature_face (const string &name) {
if (name=="linear")
};
-
string int_to_string (const unsigned int i, const unsigned int digits) {
string s;
switch (digits)
};
-
-
-
// explicit instantiations
template class FEHelper<2>;
/* $Id$ */
-
#include <base/logstream.h>
#include <lac/sparse_matrix.h>
#include <lac/vector.h>
#include <dofs/dof_constraints.h>
-
-
template <int dim>
WaveProblem<dim>::WaveProblem ()
{};
-
template <int dim>
WaveProblem<dim>::~WaveProblem ()
{};
-
template <int dim>
void WaveProblem<dim>::declare_parameters (ParameterHandler &prm)
{
};
-
template <int dim>
void WaveProblem<dim>::parse_parameters (ParameterHandler &prm)
{
};
-
template <int dim>
void WaveProblem<dim>::create_new (const unsigned int)
{
};
-
template <int dim>
void WaveProblem<dim>::run (ParameterHandler &prm)
{
// prm.print_parameters (logfile, Text);
- ////////////////////////////////
+////////////////////////////////
// Set up the time step objects
TimestepManager<dim> timestep_manager (parameters);
if (true) {
};
- ////////////////////////////////////
+////////////////////////////////////
// actually do the work (or rather:
// let the work be done)
for (unsigned int sweep=0; sweep<parameters.number_of_sweeps; ++sweep)
};
-
-
int main ()
{
// no additional output to console
// abort
return 3;
};
-
-
- return 0;
-};
+return 0;
+};
+
projects / dealii-svn.git / commitdiff