Dear Daniel,
I'm sorry to make this mistake, (and thank Jean-Paul very much). The newest
code focus on build the sparsity pattern in "setup_dof()", while other
parts (assemble_system(), assemble_rhs(), run() ) need to be left to make
the code work and maintain the structure of 3 components of the blockvector
solution.
Thank you 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 ();
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;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> system_matrix_1;
BlockSparseMatrix<double> free_energy_matrix;
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;
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
{
double phi_value = 0;
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;
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_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);
}
{
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);
}
}
template <int dim>
void Problem<dim>::assemble_system (const BlockVector<double> old_solution,
const BlockVector<double> older_solution)
{
system_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();
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::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));
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)
{
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
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));
}
}
}
}
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);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
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> > 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);
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);
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;
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];
}
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);
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);
}
}
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;
QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
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);
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::vector<double> phi_0_values(n_q_points);
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();
for (; cell!=endc; ++cell, ++cell_phi)
{
fe_values.reinit (cell);
fe_phi_values.reinit (cell_phi);
local_matrix = 0;
local_rhs = 0;
fe_phi_values.get_function_values (phi_0, phi_0_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double old_phi = phi_0_values[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 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 * (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);
}
}
}
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);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<double> phi_0_values(n_q_points);
std::vector<double> q_0_values(n_q_points);
const SolutionsPhi<dim> neumann_phi;
const SolutionsMu<dim> neumann_mu;
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();
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_q = fe_values[intervarq].value (i, q);
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);
}
}
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);
}
}
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();
FullMatrix<double> local_free_energy_matrix (dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
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);
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;
fe_values.get_function_values (solution, solution_values);
fe_values[phase].get_function_gradients(solution, solution_gradients_0);
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>::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);
}
setup_dofs ();
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);
assemble_system_1 (phi_0);
assemble_rhs_1 (phi_0, q_0);
solve_1 ();
older_solution = solution_1;
old_solution = solution_1;
assemble_system (old_solution, older_solution);
time = time_step * 2;
timestep_number = 2;
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);
}
}
}
}
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;
}
