Dear Daniel,
Thank you very much! The LittleCase0-3.cc is my work case with the same
problem as I say.
Simply describe what this code work for: to solve
<https://lh3.googleusercontent.com/-otIsCy2Aluc/W2rp2KITMnI/AAAAAAAAAFg/VhmiSjc6wR84yFJK2GHPoXHe0f-imAzsACLcBGAs/s1600/44.JPG>
from
<https://lh3.googleusercontent.com/-skgErxg4l3c/W2rp8XRaLiI/AAAAAAAAAFk/UOrruU3ZKoYUtVJpQ1TLaOT8_EszwakBQCLcBGAs/s1600/45.JPG>
where the three variables have periodic boundary conditions on the domain
[0,1]*[0,1]. So the weak form has no part of boundary integral.
Thank you very much!
Best,
Chucui
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <math.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include <deal.II/lac/sparse_direct.h>
// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include <deal.II/lac/sparse_ilu.h>
// This is C++:
#include <iostream>
#include <fstream>
#include <sstream>
// As in all programs, the namespace dealii is included:
namespace Step22
{
using namespace dealii;
template <int dim>
class Problem
{
public:
Problem (const unsigned int degree);
void run ();
private:
void setup_dofs ();
void assemble_system (const BlockVector<double> old_solution,
const BlockVector<double> older_solution);
void assemble_rhs (const double time,
const BlockVector<double> old_solution,
const BlockVector<double> older_solution);
void solve ();
void output_results_in_routine (const unsigned int timestep_number) const;
void output_results (const unsigned int refinement_cycle) const;
void output_results_phi () const;
void output_results_q () const;
void output_results_initial () const;
void refine_mesh ();
void process_solution(const double time);
double free_energy_compute (const BlockVector<double> solution);
void assemble_system_1 (const Vector<double> phi_0);
void assemble_rhs_1 (const Vector<double> phi_0, const Vector<double> q_0);
void solve_1 ();
const unsigned int degree;
Triangulation<dim> triangulation;
Triangulation<dim> triangulation_active;
FESystem<dim> fe;
FE_Q<dim> fe_phi;
FE_Q<dim> fe_mu;
FE_DGQ<dim> fe_q;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_0;
DoFHandler<dim> dof_handler_1;
DoFHandler<dim> dof_handler_2;
ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;
BlockSparsityPattern sparsity_pattern;
//BlockSparsityPattern sparsity_pattern_mu;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> system_matrix_1;
BlockSparseMatrix<double> free_energy_matrix;
//SparseMatrix<double> system_matrix_mu_0;
Vector<double> phi_0, q_0;
BlockVector<double> solution;
BlockVector<double> old_solution;
BlockVector<double> older_solution;
BlockVector<double> system_rhs;
BlockVector<double> system_rhs_1;
BlockVector<double> solution_1;
double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;
ConvergenceTable convergence_table;
};
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class SolutionsPhi : public Function<dim>
{
public:
SolutionsPhi () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double SolutionsPhi<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
//if (component == 0)
double phi_value = 0;
// phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// + std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// - (std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01 && std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01);
Tensor<1,dim> r;
r[0] = std::sqrt((p[0] - 0.3 + 0.01) * (p[0] - 0.3 + 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
r[1] = std::sqrt((p[0] - 0.7 - 0.01) * (p[0] - 0.7 - 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
if ((r[0] >= (0.2 + 0.01)) && (r[1] >= (0.2 + 0.01)))
{phi_value = -1;
}
else if (r[0] < (0.2 + 0.01))
{phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
}
else if (r[1] < (0.2 + 0.01))
{phi_value = std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
}
else
{phi_value = 0;
}
//phi_value = std::cos(numbers::PI * p[0]) * std::cos(numbers::PI * p[1]);
return phi_value;
}
template <int dim>
Tensor<1,dim> SolutionsPhi<dim>::gradient (const Point<dim> &p,
const unsigned int) const
{
Tensor<1,dim> return_value;
for (unsigned int d=0; d<dim; ++d)
{
return_value[d] = 0;}
return return_value;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class SolutionsMu : public Function<dim>
{
public:
SolutionsMu () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double SolutionsMu<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
//if (component == 0)
double phi_value = 0;
// phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// + std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// - (std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01 && std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01);
return phi_value;
}
template <int dim>
Tensor<1,dim> SolutionsMu<dim>::gradient (const Point<dim> &p,
const unsigned int) const
{
Tensor<1,dim> return_value;
for (unsigned int d=0; d<dim; ++d)
{
return_value[d] = 0;}
return return_value;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class SolutionsQ : public Function<dim>
{
public:
SolutionsQ () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double SolutionsQ<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double phi_value = 0, q_value = 0;
//phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// + std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01)*(std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)) < 0.2+0.01)
// - (std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01 && std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5))>=0.2+0.01);
Tensor<1,dim> r;
r[0] = std::sqrt((p[0] - 0.3 + 0.01) * (p[0] - 0.3 + 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
r[1] = std::sqrt((p[0] - 0.7 - 0.01) * (p[0] - 0.7 - 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
if ((r[0] >= (0.2 + 0.01)) && (r[1] >= (0.2 + 0.01)))
{phi_value = -1;
}
else if (r[0] < (0.2 + 0.01))
{phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
}
else if (r[1] < (0.2 + 0.01))
{phi_value = std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
}
else
{phi_value = 0;
}
//phi_value = std::cos(numbers::PI * p[0]) * std::cos(numbers::PI * p[1]);
q_value = (1./std::sqrt(2)) * (phi_value * phi_value - 1);
return q_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class Solutions : public Function<dim>
{
public:
Solutions () : Function<dim>(3) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double Solutions<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
if (component == 0)
return std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time());
else if (component == 1)
return std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time()) + std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time()) * std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time()) * std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time());
else
return (1 / 2) - (1 / 2) * std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time()) * std::cos(p[0]) * std::cos(p[1]) * std::cos(this->get_time());
}
template <int dim>
void
Solutions<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = Solutions<dim>::value (p, c);
}
double GFunction (const double phi1,
const double phi2)
{ const double eps = 0.01;
return (1.5 * phi1 - 0.5 * phi2) * (std::sqrt(2));
}
template <int dim>
Problem<dim>::Problem (const unsigned int degree)
:
degree (degree),
fe (FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1,
FE_DGQ<dim>(degree), 1),
fe_phi (degree),
fe_mu (degree),
fe_q (degree),
dof_handler (triangulation),
dof_handler_0 (triangulation),
dof_handler_1 (triangulation),
dof_handler_2 (triangulation),
time_step (0.001),
timestep_number (1),
D (1),
eps (0.01),
lambda (0.001)
{}
template <int dim>
void Problem<dim>::setup_dofs ()
{
//A_preconditioner.reset ();
system_matrix.clear ();
system_matrix_1.clear ();
free_energy_matrix.clear ();
dof_handler.distribute_dofs (fe);
dof_handler_0.distribute_dofs (fe_phi);
dof_handler_1.distribute_dofs (fe_mu);
dof_handler_2.distribute_dofs (fe_q);
//DoFRenumbering::Cuthill_McKee (dof_handler);
DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler_0);
DoFRenumbering::component_wise (dof_handler_1);
DoFRenumbering::component_wise (dof_handler_2);
std::vector<unsigned int> block_component (3);
block_component[0] = 0;
block_component[1] = 1;//=3-1
block_component[2] = 2;
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, component_mask);
//constraints.close ();
constraints.close ();
constraints_phi.close ();
constraints_mu.close ();
constraints_q.close ();
std::vector<types::global_dof_index> dofs_per_block (3);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
const unsigned int n_phi = dofs_per_block[0],
n_mu = dofs_per_block[1],
n_q = dofs_per_block[2];
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_phi << '+' << n_mu << '+' << n_q << ')'
<< std::endl;
{
BlockDynamicSparsityPattern dsp (3,3);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_mu, n_phi);
dsp.block(2,0).reinit (n_q, n_phi);
dsp.block(0,1).reinit (n_phi, n_mu);
dsp.block(1,1).reinit (n_mu, n_mu);
dsp.block(2,1).reinit (n_q, n_mu);
dsp.block(0,2).reinit (n_phi, n_q);
dsp.block(1,2).reinit (n_mu, n_q);
dsp.block(2,2).reinit (n_q, n_q);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
//std::ofstream out ("sparsity_pattern2.svg");
//sparsity_pattern.print_svg (out);
}
// Finally, the system matrix, solution and right hand side are created
// from the block structure as in step-20:
{
system_matrix.reinit (sparsity_pattern);
system_matrix_1.reinit (sparsity_pattern);
free_energy_matrix.reinit (sparsity_pattern);
solution.reinit (3);
solution.block(0).reinit (n_phi);
solution.block(1).reinit (n_mu);
solution.block(2).reinit (n_q);
solution.collect_sizes ();
old_solution.reinit (3);
old_solution.block(0).reinit (n_phi);
old_solution.block(1).reinit (n_mu);
old_solution.block(2).reinit (n_q);
old_solution.collect_sizes ();
older_solution.reinit (3);
older_solution.block(0).reinit (n_phi);
older_solution.block(1).reinit (n_mu);
older_solution.block(2).reinit (n_q);
older_solution.collect_sizes ();
system_rhs.reinit (3);
system_rhs.block(0).reinit (n_phi);
system_rhs.block(1).reinit (n_mu);
system_rhs.block(2).reinit (n_q);
system_rhs.collect_sizes ();
system_rhs_1.reinit (3);
system_rhs_1.block(0).reinit (n_phi);
system_rhs_1.block(1).reinit (n_mu);
system_rhs_1.block(2).reinit (n_q);
system_rhs_1.collect_sizes ();
solution_1.reinit (3);
solution_1.block(0).reinit (n_phi);
solution_1.block(1).reinit (n_mu);
solution_1.block(2).reinit (n_q);
solution_1.collect_sizes ();
phi_0.reinit (n_phi);
q_0.reinit (n_q);
}
/* {
system_matrix_mu_0.clear ();
DynamicSparsityPattern dsp_mu (n_mu, n_mu);
DoFTools::make_sparsity_pattern (dof_handler_1, dsp_mu, constraints_mu, false);
sparsity_pattern_mu.copy_from (dsp_mu);
system_matrix_mu_0.reinit (sparsity_pattern_mu);
//old_U_n.reinit (n_mu);
//U_n.reinit (n_mu);
mu_0.reinit (n_mu);
system_rhs_mu_0.reinit (n_mu);
}*/
}
template <int dim>
void Problem<dim>::assemble_system (const BlockVector<double> old_solution,
const BlockVector<double> older_solution)
{
system_matrix=0;
//std::cout << "dof 0 " << std::endl;
QGauss<dim> quadrature_formula(degree+2);
//std::cout << "dof 1 " << std::endl;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
//std::cout << "dof 2 " << std::endl;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
//std::cout << "dof 3 " << std::endl;
//const RightHandSide<dim> right_hand_side;
//std::vector<Vector<double> > rhs_values (n_q_points,
// Vector<double>(3));
std::vector<Vector<double> > older_solution_values(n_q_points, Vector<double>(3));
std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(3));
// std::vector<std::vector<Tensor<1, dim> > > older_solution_gradients(n_q_points, std::vector<Tensor<1,dim> > (3));
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar chemipt (1);
const FEValuesExtractors::Scalar intervarq (2);
//std::cout << "dof 4 " << std::endl;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
//right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
// rhs_values);
fe_values.get_function_values (old_solution, old_solution_values);
fe_values.get_function_values (older_solution, older_solution_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double older_phi = older_solution_values[q](0);
const double old_phi = old_solution_values[q](0);
const Tensor<1,dim> grad_psi_i_phi = fe_values[phase].gradient (i, q);
const Tensor<1,dim> grad_psi_i_mu = fe_values[chemipt].gradient (i, q);
const double psi_i_phi = fe_values[phase].value (i, q);
const double psi_i_mu = fe_values[chemipt].value (i, q);
const double psi_i_q = fe_values[intervarq].value (i, q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> grad_psi_j_phi = fe_values[phase].gradient (j, q);
const Tensor<1,dim> grad_psi_j_mu = fe_values[chemipt].gradient (j, q);
const double psi_j_phi = fe_values[phase].value (j, q);
const double psi_j_mu = fe_values[chemipt].value (j, q);
const double psi_j_q = fe_values[intervarq].value (j, q);
local_matrix(i,j) += (psi_i_phi * psi_j_phi
+ time_step * 0.5 * D * grad_psi_i_phi * grad_psi_j_mu
+ psi_i_mu * psi_j_mu
- grad_psi_i_mu * grad_psi_j_phi * eps * lambda
- (lambda / eps) * GFunction(old_phi, older_phi) * psi_i_mu * psi_j_q
+ psi_i_q * psi_j_q
- GFunction(old_phi, older_phi) * psi_i_q * psi_j_phi
) * fe_values.JxW(q);
}
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
local_matrix(i,j));
}
}
}
//std::cout << " Computing preconditioner..." << std::endl << std::flush;
}
template <int dim>
void Problem<dim>::assemble_rhs (const double time,
const BlockVector<double> old_solution,
const BlockVector<double> older_solution)
{
system_rhs=0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
//std::cout << "dof 9-1 " << std::endl;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//std::cout << "dof 9-2 " << std::endl;
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
//std::cout << "dof 9-3 " << std::endl;
std::vector<Vector<double> > older_solution_values(n_q_points, Vector<double>(3));
std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(3));
std::vector<Tensor<1, dim> > old_solution_gradients_0(n_q_points);
std::vector<Tensor<1, dim> > old_solution_gradients_1(n_q_points);
//std::vector<std::vector<Tensor<1, dim> > > old_solution_gradients(n_q_points, std::vector<Tensor<1,dim> >(3));
//std::cout << "dof 9-5 " << std::endl;
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar chemipt (1);
const FEValuesExtractors::Scalar intervarq (2);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
local_rhs = 0;
fe_values.reinit (cell);
fe_values.get_function_values (old_solution, old_solution_values);
fe_values.get_function_values (older_solution, older_solution_values);
fe_values[phase].get_function_gradients(old_solution, old_solution_gradients_0);
fe_values[chemipt].get_function_gradients(old_solution, old_solution_gradients_1);
//right_hand_side.vector_value_list(fe_values.get_quadrature_points(),rhs_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double older_phi = older_solution_values[q](0);
const double old_phi = old_solution_values[q](0);
const double old_mu = old_solution_values[q](1);
const double old_q = old_solution_values[q](2);
Tensor<1,dim> grad_old_phi;
Tensor<1,dim> grad_old_mu;
//std::vector<Vector<double> > grad_old(n_q_points, Vector<double>(3));
for (unsigned int d=0; d<dim; ++d)
{
grad_old_phi[d] = old_solution_gradients_0[q][d];
grad_old_mu[d] = old_solution_gradients_1[q][d];
}
//const double new_fright = rhs_values[q](0);
//const double new_fright = 0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> grad_psi_i_phi = fe_values[phase].gradient (i, q);
const Tensor<1,dim> grad_psi_i_mu = fe_values[chemipt].gradient (i, q);
const double psi_i_phi = fe_values[phase].value (i, q);
const double psi_i_mu = fe_values[chemipt].value (i, q);
const double psi_i_q = fe_values[intervarq].value (i, q);
//const unsigned int component_i = fe.system_to_component_index(i).first;
local_rhs(i) += (//psi_i_phi * new_fright * time_step
+ psi_i_phi * old_phi
- 0.5 * time_step * D * grad_psi_i_phi * grad_old_mu
- psi_i_mu * old_mu
+ grad_psi_i_mu * grad_old_phi * eps * lambda
+ (lambda / eps) * GFunction(old_phi, older_phi) * psi_i_mu * old_q
+ psi_i_q * old_q
- GFunction(old_phi, older_phi) * psi_i_q * old_phi)
* fe_values.JxW(q);
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
//std::cout << "dof 9-6 " << std::endl;
}
template <int dim>
void Problem<dim>::solve ()
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult (solution, system_rhs);
}
template <int dim>
void Problem<dim>::assemble_system_1 (const Vector<double> phi_0)
{
system_matrix_1=0;
//std::cout << "dof 0 " << std::endl;
QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
//std::cout << "dof 1 " << std::endl;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
FEValues<dim> fe_phi_values (fe_phi, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
//std::cout << "dof 2 " << std::endl;
// FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
// update_values | update_normal_vectors | update_gradients |
// update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
//const unsigned int n_face_q_points = face_quadrature_formula.size();
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
//std::cout << "dof 3 " << std::endl;
std::vector<double> phi_0_values(n_q_points);
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar chemipt (1);
const FEValuesExtractors::Scalar intervarq (2);
//std::cout << "dof 4 " << std::endl;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
cell_phi = dof_handler_0.begin_active();
//std::cout << "dof 12-0-16 " << std::endl;
for (; cell!=endc; ++cell, ++cell_phi)
{
fe_values.reinit (cell);
fe_phi_values.reinit (cell_phi);
local_matrix = 0;
local_rhs = 0;
//right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
// rhs_values);
fe_phi_values.get_function_values (phi_0, phi_0_values);
//fe_values.get_function_values (older_solution, older_solution_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double old_phi = phi_0_values[q];
//std::cout << "dof 12-0 " << old_phi << std::endl;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> grad_psi_i_phi = fe_values[phase].gradient (i, q);
const Tensor<1,dim> grad_psi_i_mu = fe_values[chemipt].gradient (i, q);
const double psi_i_phi = fe_values[phase].value (i, q);
const double psi_i_mu = fe_values[chemipt].value (i, q);
const double psi_i_q = fe_values[intervarq].value (i, q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> grad_psi_j_phi = fe_values[phase].gradient (j, q);
const Tensor<1,dim> grad_psi_j_mu = fe_values[chemipt].gradient (j, q);
const double psi_j_phi = fe_values[phase].value (j, q);
const double psi_j_mu = fe_values[chemipt].value (j, q);
const double psi_j_q = fe_values[intervarq].value (j, q);
//std::cout << "dof 12-0 " << "|"<< i<< "|"<< j<< "|"<< psi_i_phi * psi_j_phi * (1./time_step) << std::endl;
//std::cout << "dof 12-1 " << "|"<< i<< "|"<< j<< "|"<< +grad_psi_i_phi[0] * grad_psi_j_mu[0] + grad_psi_i_phi[1] * grad_psi_j_mu[1] << std::endl;
//std::cout << "dof 12-2 " << "|"<< i<< "|"<< j<< "|"<< +psi_i_mu * psi_j_mu << std::endl;
//std::cout << "dof 12-3 " << "|"<< i<< "|"<< j<< "|"<< -(grad_psi_i_mu[0] * grad_psi_j_phi[0] + grad_psi_i_mu[1] * grad_psi_j_phi[1]) * eps * lambda << std::endl;
//std::cout << "dof 12-4 " << "|"<< i<< "|"<< j<< "|"<< -(lambda / eps) * GFunction(old_phi, old_phi) * psi_i_mu * psi_j_q << std::endl;
//std::cout << "dof 12-5 " << "|"<<i<< "|"<< j<< "|"<< +psi_i_q * psi_j_q << std::endl;
//std::cout << "dof 12-6 " << "|"<< i<< "|"<< j<< "|"<< -GFunction(old_phi, old_phi) * psi_i_q * psi_j_phi << std::endl;
local_matrix(i,j) += (psi_i_phi * psi_j_phi * (1./time_step )
+ (grad_psi_i_phi * grad_psi_j_mu)
+ psi_i_mu * psi_j_mu
- grad_psi_i_mu * grad_psi_j_phi * eps * lambda
- (lambda / eps) * GFunction(old_phi, old_phi) * psi_i_mu * psi_j_q
+ psi_i_q * psi_j_q
- GFunction(old_phi, old_phi) * psi_i_q * psi_j_phi
)
* fe_values.JxW(q);
//std::cout << "dof 12-6 " << "|"<< i<< "|"<< j<< "|"<< local_matrix(i,j) << std::endl;
}
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
system_matrix_1.add (local_dof_indices[i],
local_dof_indices[j],
local_matrix(i,j));
}
}
}
}
template <int dim>
void Problem<dim>::assemble_rhs_1 (const Vector<double> phi_0,
const Vector<double> q_0)
{
system_rhs_1=0;
QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEValues<dim> fe_phi_values (fe_phi, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEValues<dim> fe_q_values (fe_q, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
//FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
// update_values | update_normal_vectors | update_gradients |
// update_quadrature_points | update_JxW_values);
//const unsigned int n_face_q_points = face_quadrature_formula.size();
//std::cout << "dof 9-1 " << std::endl;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//std::cout << "dof 9-2 " << std::endl;
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
//std::cout << "dof 9-3 " << std::endl;
std::vector<double> phi_0_values(n_q_points);//, Vector<double>(3));
std::vector<double> q_0_values(n_q_points);//, Vector<double>(3));
//std::vector<Tensor<1, dim> > old_solution_gradients_0(n_q_points);
//std::vector<Tensor<1, dim> > old_solution_gradients_1(n_q_points);
//std::vector<std::vector<Tensor<1, dim> > > old_solution_gradients(n_q_points, std::vector<Tensor<1,dim> >(3));
const SolutionsPhi<dim> neumann_phi;
const SolutionsMu<dim> neumann_mu;
//std::cout << "dof 9-5 " << std::endl;
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar chemipt (1);
const FEValuesExtractors::Scalar intervarq (2);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
cell_phi = dof_handler_0.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_q = dof_handler_2.begin_active();
//std::cout << "dof 12-0-16 " << std::endl;
for (; cell!=endc; ++cell, ++cell_phi, ++cell_q)
{
local_rhs = 0;
fe_values.reinit (cell);
fe_phi_values.reinit (cell_phi);
fe_q_values.reinit (cell_q);
fe_phi_values.get_function_values (phi_0, phi_0_values);
fe_q_values.get_function_values (q_0, q_0_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double old_phi = phi_0_values[q];
const double old_q = q_0_values[q];
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double psi_i_phi = fe_values[phase].value (i, q);
//const double psi_i_mu = fe_values[chemipt].value (i, q);
const double psi_i_q = fe_values[intervarq].value (i, q);
//const unsigned int component_i = fe.system_to_component_index(i).first;
local_rhs(i) += (psi_i_phi * old_phi * (1./time_step)
+ psi_i_q * old_q
- GFunction(old_phi, old_phi) * psi_i_q * old_phi)
* fe_values.JxW(q);
//std::cout << "dof 12-6 " << "|"<< i<< "|"<< local_rhs(i) << std::endl;
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs_1(local_dof_indices[i]) += local_rhs(i);
}
//std::cout << "dof 9-6 " << std::endl;
}
template <int dim>
void Problem<dim>::solve_1 ()
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix_1);
A_direct.vmult (solution_1, system_rhs_1);
}
template <int dim>
double Problem<dim>::free_energy_compute (const BlockVector<double> solution)
{
free_energy_matrix=0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//std::cout << "n_q_points " << n_q_points << std::endl;
FullMatrix<double> local_free_energy_matrix (dofs_per_cell, dofs_per_cell);
//Vector<double> local_rhs (dofs_per_cell);
//std::cout << "dofs_per_cell " << dofs_per_cell << std::endl;
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
///std::vector<Vector<double> > older_solution_values(n_q_points, Vector<double>(3));
std::vector<Vector<double> > solution_values(n_q_points, Vector<double>(3));
std::vector<Tensor<1, dim> > solution_gradients_0(n_q_points);
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar chemipt (1);
const FEValuesExtractors::Scalar intervarq (2);
//std::cout << "dof 4 " << std::endl;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
local_free_energy_matrix = 0;
//right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
// rhs_values);
fe_values.get_function_values (solution, solution_values);
fe_values[phase].get_function_gradients(solution, solution_gradients_0);
//fe_values.get_function_values (older_solution, older_solution_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> grad_psi_i_phi = fe_values[phase].gradient (i, q);
const double psi_i_q = fe_values[intervarq].value (i, q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> grad_psi_j_phi = fe_values[phase].gradient (j, q);
const double psi_j_q = fe_values[intervarq].value (j, q);
local_free_energy_matrix(i,j) += (0.5 * lambda * eps * grad_psi_i_phi * grad_psi_j_phi
+ psi_i_q * psi_j_q * 0.5 * (lambda / eps))
* fe_values.JxW(q);
}
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
free_energy_matrix.add (local_dof_indices[i],
local_dof_indices[j],
local_free_energy_matrix(i,j));
}
}
}
double free_energy = 0;
free_energy = free_energy_matrix.matrix_norm_square(solution);
//std::cout << "dof 11-4-8 "<< free_energy << std::endl;
return free_energy;
}
template <int dim>
void Problem<dim>::output_results_initial () const
{
std::vector<std::string> solution_names;
solution_names.push_back ("phi");
solution_names.push_back ("mu");
solution_names.push_back ("q");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution_1, solution_names);
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "solution_1"
//<< Utilities::int_to_string(timestep_number)
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Problem<dim>::output_results_in_routine (const unsigned int timestep_number) const
{
std::vector<std::string> solution_names;
solution_names.push_back ("phi");
solution_names.push_back ("mu");
solution_names.push_back ("q");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (old_solution, solution_names);
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "solution-"
<< Utilities::int_to_string(timestep_number)
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Problem<dim>::output_results_phi () const
{
std::vector<std::string> solution_names;
solution_names.push_back ("phi");
//solution_names.push_back ("q");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler_0);
data_out.add_data_vector (phi_0, solution_names);
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "phi-0"
//<< Utilities::int_to_string()
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Problem<dim>::output_results_q () const
{
std::vector<std::string> solution_names;
//solution_names.push_back ("phi");
solution_names.push_back ("q");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler_2);
data_out.add_data_vector (q_0, solution_names);
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "q-0"
//<< Utilities::int_to_string()
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Problem<dim>::output_results (const unsigned int refinement_cycle) const
{
//if (timestep_number % 5 != 0)
// return;
std::vector<std::string> solution_names;
solution_names.push_back ("phi");
solution_names.push_back ("mu");
solution_names.push_back ("q");
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names);
data_out.build_patches (degree+3);
std::ostringstream filename;
filename << "solution-"
<< Utilities::int_to_string(refinement_cycle,3)
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Problem<dim>::process_solution(const double time)
{
Solutions<dim> exact_solution;
SolutionsPhi<dim> exact_solution_phi;
exact_solution_phi.set_time (time);
Vector<float> difference_per_cell(triangulation.n_active_cells());
VectorTools::integrate_difference(dof_handler_0,
solution.block(0),
exact_solution_phi,
difference_per_cell,
QGauss<dim>(degree+2),
VectorTools::L2_norm);
const double L2_error = VectorTools::compute_global_error(
triangulation, difference_per_cell, VectorTools::L2_norm);
VectorTools::integrate_difference(dof_handler_0,
solution.block(0),
exact_solution_phi,
difference_per_cell,
QGauss<dim>(degree+2),
VectorTools::H1_seminorm);
const double H1_error = VectorTools::compute_global_error(
triangulation, difference_per_cell, VectorTools::H1_seminorm);
const QTrapez<1> q_trapez;
const QIterated<dim> q_iterated(q_trapez, 5);
VectorTools::integrate_difference(dof_handler_0,
solution.block(0),
exact_solution_phi,
difference_per_cell,
q_iterated,
VectorTools::Linfty_norm);
const double Linfty_error = VectorTools::compute_global_error(
triangulation, difference_per_cell, VectorTools::Linfty_norm);
const unsigned int n_active_cells = triangulation.n_active_cells();
const unsigned int n_dofs = dof_handler.n_dofs();
//std::cout << "Cycle " << cycle << ':' << std::endl
// << " Number of active cells: " << n_active_cells
// << std::endl
// << " Number of degrees of freedom: " << n_dofs << std::endl;
//convergence_table.add_value("cycle", cycle);
convergence_table.add_value("cells", n_active_cells);
convergence_table.add_value("dofs", n_dofs);
convergence_table.add_value("L2", L2_error);
convergence_table.add_value("H1", H1_error);
convergence_table.add_value("Linfty", Linfty_error);
}
template <int dim>
void Problem<dim>::run ()
{
const unsigned int n_cycles = 1;
double free_energy = 0;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation_active, 0, 1);
triangulation_active.refine_global (7);
GridGenerator::flatten_triangulation (triangulation_active, triangulation);
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary())
{
if (std::fabs(cell->face(f)->center()(0) - (1)) < 1e-12)
cell->face(f)->set_boundary_id (1);
else if(std::fabs(cell->face(f)->center()(1) - (0)) < 1e-12)
cell->face(f)->set_boundary_id (2);
else if(std::fabs(cell->face(f)->center()(1) - (1)) < 1e-12)
cell->face(f)->set_boundary_id (3);
else if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
{
cell->face(f)->set_boundary_id (0);
}
}
}
}
}
else
{
//triangulation.refine_global (1);
}
std::cout << "before_start" << cycle << std::endl;
setup_dofs ();
// std::cout << "start" << cycle << std::endl;
//std::cout << "phi_size" << dof_handler_0.n_dofs() << phi_0.size() << std::endl;
// std::cout << "q_size" << dof_handler_2.n_dofs() << q_0.size() << std::endl;
VectorTools::project (dof_handler_0, constraints_phi, QGauss<dim>(degree+2),
SolutionsPhi<dim>(),
phi_0);
VectorTools::project (dof_handler_2, constraints_q, QGauss<dim>(degree+2),
SolutionsQ<dim>(),
q_0);
output_results_phi ();
output_results_q ();
// std::cout << "initial" << cycle << std::endl;
assemble_system_1 (phi_0);
// std::cout << "system" << cycle << std::endl;
assemble_rhs_1 (phi_0, q_0);
solve_1 ();
// std::cout << "solve" << cycle << std::endl;
older_solution = solution_1;
old_solution = solution_1;
output_results_initial ();
assemble_system (old_solution, older_solution);
time = time_step * 2;
timestep_number = 2;
//std::cout << "Time step start " << timestep_number << " at t=" << time << std::endl;
for (; time<=300.01; time+=time_step, ++timestep_number) //time<=/ 1./16
{
assemble_rhs (time - 0.5 * time_step, old_solution, older_solution);
solve ();
free_energy = free_energy_compute (solution);
older_solution = old_solution;
old_solution = solution;
std::cout << "Time step " << timestep_number << " at t=" << time << "free energy is " << free_energy << std::endl;
assemble_system (old_solution, older_solution);
output_results_in_routine (timestep_number);
}
////compute_errors ();
output_results(cycle);
// process_solution(time - 1.0 * time_step);
}
/* convergence_table.set_precision("L2", 3);
convergence_table.set_precision("H1", 3);
convergence_table.set_precision("Linfty", 3);
convergence_table.set_scientific("L2", true);
convergence_table.set_scientific("H1", true);
convergence_table.set_scientific("Linfty", true);
convergence_table.set_tex_caption("cells", "\\# cells");
convergence_table.set_tex_caption("dofs", "\\# dofs");
convergence_table.set_tex_caption("L2", "@f$L^2@f$-error");
convergence_table.set_tex_caption("H1", "@f$H^1@f$-error");
convergence_table.set_tex_caption("Linfty", "@f$L^\\infty@f$-error");
convergence_table.set_tex_format("cells", "r");
convergence_table.set_tex_format("dofs", "r");
std::cout << std::endl;
//error_filename += ".tex";
//std::ofstream error_table_file(error_filename.c_str());
convergence_table.write_text(std::cout);
std::string error_filename = "error";
error_filename += ".tex";*/
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
Problem<2> flow_problem(1);
flow_problem.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
