Dear Jean-Paul,
Thank you very much for your patience, when I print some information like
this:
template <int dim>
void Problem<dim>::run ()
{
[...]
setup_dofs ();
std::cout << "test-1 "<< std::endl;
VectorTools::project (dof_handler_0, constraints_phi,
QGauss<dim>(degree+2),
SolutionsPhi<dim>(),
phi_0);
std::cout << "test-2 " << std::endl;
VectorTools::project (dof_handler_2, constraints_q,
QGauss<dim>(degree+2),
SolutionsQ<dim>(),
q_0);
std::cout << "test-3 " << std::endl;
assemble_system_1 (phi_0);
std::cout << "test-4" << std::endl;
assemble_rhs_1 (phi_0, q_0);
std::cout << "test-5" << std::endl;
and the output is
Linking CXX executable Case0-3
[ 50%] Built target Case0-3
[100%] Run with Debug configuration
Number of active cells: 16384
Number of degrees of freedom: 98818 (16641+16641+65536)
test-1
test-2
test-3
--------------------------------------------------------
An error occurred in line <1764> of file
</home/yucan/boost/deal.ii-8.5.1/deal.II/include/deal.II/lac/sparse_matrix.h>
in function
void
dealii::SparseMatrix<number>::add(dealii::SparseMatrix<number>::size_type,
dealii::SparseMatrix<number>::size_type, number) [with number = double;
dealii::SparseMatrix<number>::size_type = unsigned int]
The violated condition was:
(index != SparsityPattern::invalid_entry) || (value == number())
Additional information:
You are trying to access the matrix entry with index <0,1>, but this
entry does not exist in the sparsity pattern of this matrix.
So it can't run when arrive at "assemble_system_1", so "setup-dof()",
"assemble_system_1()", "run()" are left—— this is LittleCase0-3.cc
if only "setup_dof()", "run()" are left——this is LittleCase0-4.cc
Hope this is a minimal example. 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_1 (const Vector<double> phi_0);
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;
}
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_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>::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 ();
std::cout << "test-1 "<< std::endl;
VectorTools::project (dof_handler_0, constraints_phi, QGauss<dim>(degree+2),
SolutionsPhi<dim>(),
phi_0);
std::cout << "test-2 " << std::endl;
std::cout << "test-3 " << std::endl;
assemble_system_1 (phi_0);
}
}
}
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;
}
#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_1 (const Vector<double> phi_0);
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>
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>::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 ();
}
}
}
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;
}
