The attachment are 2 code files：LittleCase0-4.cc (for test),
LittleCase0-5.cc (RT element).
在 2018年8月26日星期日 UTC+8下午8:32:40，[email protected]写道：
>
> Hi, All,
>
> I have a question of behavior of FESystem with 2 vectors and 2 scalars, my
> original problem is to solve:
>
> [image: 56.JPG]
> for \phi, \sigma, \mu, \gamma, I choose FESystem to describe the
> information of them
> FESystem<dim> fe;
> [...]
> fe (FE_DGQ<dim>(degree), 1,
> FE_RaviartThomas<dim>(degree), 1,
> FE_DGQ<dim>(degree), 1,
> FE_RaviartThomas<dim>(degree), 1),
>
> The test code is LittleCase0-5.cc, and after make run, the error
> information is:
> Linking CXX executable Case0-5
> [ 50%] Built target Case0-5
> [100%] Run with Debug configuration
> before_start0
> dof 9-1
> dof 9-2 6|112
> dof 9-3
> dof 9-4
> dof 9-5
>
> --------------------------------------------------------
> An error occurred in line <1746> of file
> </home/yucan/boost/deal.ii-8.5.1/src/source/dofs/dof_tools.cc> in function
> void dealii::DoFTools::count_dofs_per_block(const DoFHandlerType&,
> std::vector<unsigned int>&, const std::vector<unsigned int>&) [with
> DoFHandlerType = dealii::DoFHandler<2>]
> The violated condition was:
> target_block.size()==fe.n_blocks()
> Additional information:
> Dimension 6 not equal to 4.
>
>
> to test the behavior of FESystem with 2 vectors and 2 scalars, I choose
> another FESystem:
> FESystem<dim> fe;
> [...]
> fe (FE_Q<dim>(degree), 1,
> FE_Q<dim>(degree), dim,
> FE_Q<dim>(degree), 1,
> FE_Q<dim>(degree), dim),
>
> the code is LittleCase0-4.cc, and the result after "make run" is:
> Linking CXX executable Case0-4
> [ 50%] Built target Case0-4
> [100%] Run with Debug configuration
> before_start0
> dof 9-1
> dof 9-2 6|54
> dof 9-3
> dof 9-4
> dof 9-5
> dof 9-6
> dof 9-7
> Number of active cells: 4
> Number of degrees of freedom: 54 (9+27|9+9)
> dof 9-8
> dof 9-9
> dof 9-10
> dof 9-11
> dof 9-12
> dof 9-13
> setup_dofs0
> [100%] Built target run
>
>
> My question is:
>
> 1. Can FESystem countain 2 vectors and 2 scalars? If so, how to make
> "dof_handler", "dofs_per_block", and "block_component (target_block)"
> correctly matching? As above, I make some faults. For example, in
> LittleCase0-4.cc, I use this defination for "block_component"
> (target_block):
> std::vector<unsigned int> block_component (2*dim+2,1);
> block_component[0] = 0;
> block_component[1] = 1;
> block_component[dim+1] = 2;//=3-1
> block_component[dim+2] = 3;
> and te result is
> Number of degrees of freedom: 54 (9+27|9+9)
> but it should be: " 54 (9+18 | 9+18) ", I don't know how to control it as
> I cannot set "std::vector<unsigned int> block_component (2*dim+2,1,3);",
> where "1" and "3" means that the vector block should be in the first and
> third block in the whole.
>
> 2. Can I use RT element by this way (in Block)? If so, how to correct the
> LittleCase0-5.cc? If not, how to implement the original problem? Maybe use
> 2 FESystem fe_1, fe_2, both of which contain only 1 vector (RT element) and
> 1 scalar (DG element)?
>
>
> Thanks in advance!
>
> 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/fe/fe_raviart_thomas.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>
#include<stdlib.h>
#include "stdio.h"
#include<time.h>
// This is C++:
#include <iostream>
#include <fstream>
#include <sstream>
#define random(x) (rand()%x)
#define MAX 2
// As in all programs, the namespace dealii is included:
namespace Step22
{
using namespace dealii;
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int degree);
void run ();
private:
void setup_dofs ();
void refine_mesh ();
const unsigned int degree;
Triangulation<dim> triangulation_active;
Triangulation<dim> triangulation;
FESystem<dim> fe;
FE_Q<dim> fe_phi;
FE_Q<dim> fe_mu;
FE_Q<dim> fe_q;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_phi;
DoFHandler<dim> dof_handler_mu;
DoFHandler<dim> dof_handler_q;
ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;
BlockSparsityPattern sparsity_pattern;
SparsityPattern sparsity_pattern_q;
BlockSparseMatrix<double> system_matrix;
double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;
ConvergenceTable convergence_table;
};
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
:
degree (degree),
fe (FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), dim,
FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), dim),
fe_phi (degree),
fe_mu (degree),
fe_q (degree),
dof_handler (triangulation),
dof_handler_phi (triangulation),
dof_handler_mu (triangulation),
dof_handler_q (triangulation),
time_step (0.001),
timestep_number (1),
D (1),
eps (0.01),
lambda (0.001)
{}
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
//A_preconditioner.reset ();
system_matrix.clear ();
dof_handler.distribute_dofs (fe);
DoFRenumbering::Cuthill_McKee (dof_handler);
std::cout << "dof 9-1 " << std::endl;
std::vector<unsigned int> block_component (2*dim+2,1);
block_component[0] = 0;
block_component[1] = 1;
block_component[dim+1] = 2;//=3-1
block_component[dim+2] = 3;
//block_component[2] = 2;//=3-1
//block_component[3] = 3;
std::cout << "dof 9-2 " << block_component.size() << "|" << dof_handler.n_dofs() << std::endl;
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
std::cout << "dof 9-3 " << std::endl;
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, component_mask);
std::cout << "dof 9-4 " << std::endl;
constraints.close ();
std::vector<types::global_dof_index> dofs_per_block (4);
std::cout << "dof 9-5 " << std::endl;
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
std::cout << "dof 9-6 " << std::endl;
const unsigned int n_phi = dofs_per_block[0],
n_gamma = dofs_per_block[1],
n_mu = dofs_per_block[2],
n_sigma = dofs_per_block[3];
std::cout << "dof 9-7 " << std::endl;
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_phi << "+" << n_gamma << "|" << n_mu << "+" << n_sigma << ")"
<< std::endl;
std::cout << "dof 9-8 " << std::endl;
dof_handler_phi.distribute_dofs (fe_phi);
std::cout << "dof 9-9 " << std::endl;
DynamicSparsityPattern dsp_phi (n_phi, n_phi);
std::cout << "dof 9-10 " << std::endl;
constraints_phi.clear ();
std::cout << "dof 9-11 " << std::endl;
DoFTools::make_sparsity_pattern (dof_handler_phi, dsp_phi, constraints_phi, false);
std::cout << "dof 9-12 " << std::endl;
DoFRenumbering::component_wise (dof_handler_phi);
std::cout << "dof 9-13 " << std::endl;
constraints_phi.close ();
}
template <int dim>
void StokesProblem<dim>::run ()
{
const unsigned int n_cycles = 1;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation_active, 0, 1.0);
triangulation_active.refine_global (1);
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.0)) < 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.0)) < 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 << "setup_dofs" << cycle << std::endl;
}
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
StokesProblem<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/fe/fe_raviart_thomas.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>
#include<stdlib.h>
#include "stdio.h"
#include<time.h>
// This is C++:
#include <iostream>
#include <fstream>
#include <sstream>
#define random(x) (rand()%x)
#define MAX 2
// As in all programs, the namespace dealii is included:
namespace Step22
{
using namespace dealii;
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int degree);
void run ();
private:
void setup_dofs ();
const unsigned int degree;
Triangulation<dim> triangulation_active;
Triangulation<dim> triangulation;
FESystem<dim> fe;
FE_DGQ<dim> fe_phi;
FE_DGQ<dim> fe_mu;
FE_DGQ<dim> fe_q;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_phi;
DoFHandler<dim> dof_handler_mu;
DoFHandler<dim> dof_handler_q;
ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;
BlockSparsityPattern sparsity_pattern;
SparsityPattern sparsity_pattern_q;
BlockSparseMatrix<double> system_matrix;
double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;
ConvergenceTable convergence_table;
};
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
:
degree (degree),
fe (FE_DGQ<dim>(degree), 1,
FE_RaviartThomas<dim>(degree), 1,
FE_DGQ<dim>(degree), 1,
FE_RaviartThomas<dim>(degree), 1),
fe_phi (degree),
fe_mu (degree),
fe_q (degree),
dof_handler (triangulation),
dof_handler_phi (triangulation),
dof_handler_mu (triangulation),
dof_handler_q (triangulation),
time_step (0.001),
timestep_number (1),
D (1),
eps (0.01),
lambda (0.001)
{}
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
//A_preconditioner.reset ();
system_matrix.clear ();
dof_handler.distribute_dofs (fe);
DoFRenumbering::Cuthill_McKee (dof_handler);
std::cout << "dof 9-1 " << std::endl;
std::vector<unsigned int> block_component (2*dim+2,1);
block_component[0] = 0;
block_component[1] = 1;
block_component[dim+1] = 2;//
block_component[dim+2] = 3;
//block_component[2] = 2;//
//block_component[3] = 3;
std::cout << "dof 9-2 " << block_component.size() << "|" << dof_handler.n_dofs() << std::endl;
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
std::cout << "dof 9-3 " << std::endl;
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, component_mask);
std::cout << "dof 9-4 " << std::endl;
constraints.close ();
std::vector<types::global_dof_index> dofs_per_block (4);
std::cout << "dof 9-5 " << std::endl;
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
std::cout << "dof 9-6 " << std::endl;
const unsigned int n_phi = dofs_per_block[0],
n_gamma = dofs_per_block[1],
n_mu = dofs_per_block[2],
n_sigma = dofs_per_block[3];
std::cout << "dof 9-7 " << std::endl;
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_phi << "+" << n_gamma << "|" << n_mu << "+" << n_sigma << ")"
<< std::endl;
std::cout << "dof 9-8 " << std::endl;
dof_handler_phi.distribute_dofs (fe_phi);
std::cout << "dof 9-9 " << std::endl;
DynamicSparsityPattern dsp_phi (n_phi, n_phi);
std::cout << "dof 9-10 " << std::endl;
constraints_phi.clear ();
std::cout << "dof 9-11 " << std::endl;
DoFTools::make_sparsity_pattern (dof_handler_phi, dsp_phi, constraints_phi, false);
std::cout << "dof 9-12 " << std::endl;
DoFRenumbering::component_wise (dof_handler_phi);
std::cout << "dof 9-13 " << std::endl;
constraints_phi.close ();
}
template <int dim>
void StokesProblem<dim>::run ()
{
const unsigned int n_cycles = 1;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation_active, 0, 1.0);
triangulation_active.refine_global (1);
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.0)) < 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.0)) < 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 << "setup_dofs" << cycle << std::endl;
}
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
StokesProblem<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;
}
