#include <grid/tria.h>
#include <grid/grid_generator.h>
+#include <grid/intergrid_map.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
const unsigned int degree;
Triangulation<dim> triangulation;
+
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
+ FE_Q<dim> scalar_velocity_fe;
+ DoFHandler<dim> scalar_velocity_dof_handler;
+
ConstraintMatrix hanging_node_constraints;
+ ConstraintMatrix scalar_velocity_hanging_node_constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
+ SparsityPattern scalar_velocity_sparsity_pattern;
+ SparseMatrix<double> scalar_velocity_system_matrix;
+ std::vector<unsigned int> scalar_to_vector_velocities[dim];
+
+
double time_step;
unsigned int timestep_number;
FE_Q<dim>(degree), 1,
FE_DGQ<dim>(degree-1), 1),
dof_handler (triangulation),
+ scalar_velocity_fe (degree+1),
+ scalar_velocity_dof_handler (triangulation),
time_step (0),
rebuild_matrices (true),
rebuild_preconditioner (true)
<< std::endl
<< std::endl;
+
+ // also initialize the DoFHandler
+ // for the scalar version of the
+ // Laplace equation on the velocity
+ // variables
+ scalar_velocity_dof_handler.distribute_dofs (scalar_velocity_fe);
+ DoFRenumbering::Cuthill_McKee (scalar_velocity_dof_handler);
+ scalar_velocity_hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (scalar_velocity_dof_handler,
+ scalar_velocity_hanging_node_constraints);
+ scalar_velocity_hanging_node_constraints.close ();
+
// The next step is to
// create the sparsity
// pattern for the system matrix
// the sparsity pattern.
if (setup_matrices == true)
{
- system_matrix.clear ();
+ {
+ system_matrix.clear ();
- BlockCompressedSetSparsityPattern csp (3,3);
+ BlockCompressedSetSparsityPattern csp (3,3);
- csp.block(0,0).reinit (n_u, n_u);
- csp.block(0,1).reinit (n_u, n_p);
- csp.block(0,2).reinit (n_u, n_T);
- csp.block(1,0).reinit (n_p, n_u);
- csp.block(1,1).reinit (n_p, n_p);
- csp.block(1,2).reinit (n_p, n_T);
- csp.block(2,0).reinit (n_T, n_u);
- csp.block(2,1).reinit (n_T, n_p);
- csp.block(2,2).reinit (n_T, n_T);
+ csp.block(0,0).reinit (n_u, n_u);
+ csp.block(0,1).reinit (n_u, n_p);
+ csp.block(0,2).reinit (n_u, n_T);
+ csp.block(1,0).reinit (n_p, n_u);
+ csp.block(1,1).reinit (n_p, n_p);
+ csp.block(1,2).reinit (n_p, n_T);
+ csp.block(2,0).reinit (n_T, n_u);
+ csp.block(2,1).reinit (n_T, n_p);
+ csp.block(2,2).reinit (n_T, n_T);
+
+ csp.collect_sizes ();
- csp.collect_sizes ();
+ DoFTools::make_sparsity_pattern (dof_handler, csp);
+ hanging_node_constraints.condense (csp);
+ sparsity_pattern.copy_from (csp);
+
+ system_matrix.reinit (sparsity_pattern);
+ }
- DoFTools::make_sparsity_pattern (dof_handler, csp);
- hanging_node_constraints.condense (csp);
- sparsity_pattern.copy_from (csp);
+ {
+ scalar_velocity_system_matrix.clear ();
- system_matrix.reinit (sparsity_pattern);
+ CompressedSetSparsityPattern csp (scalar_velocity_dof_handler.n_dofs(),
+ scalar_velocity_dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (scalar_velocity_dof_handler, csp);
+ scalar_velocity_hanging_node_constraints.condense (csp);
+ scalar_velocity_sparsity_pattern.copy_from (csp);
+
+ scalar_velocity_system_matrix.reinit (scalar_velocity_sparsity_pattern);
+ }
}
+
+ {
+ InterGridMap<DoFHandler<dim> > scalar_to_vector;
+ scalar_to_vector.make_mapping (scalar_velocity_dof_handler,
+ dof_handler);
+ std::vector<unsigned int>
+ cell_dof_indices (fe.dofs_per_cell),
+ scalar_velocity_cell_dof_indices (scalar_velocity_fe.dofs_per_cell);
+
+ for (unsigned int component=0; component<dim; ++component)
+ {
+ scalar_to_vector_velocities[component].resize (scalar_velocity_dof_handler.n_dofs());
+ std::fill (scalar_to_vector_velocities[component].begin(),
+ scalar_to_vector_velocities[component].end(),
+ numbers::invalid_unsigned_int);
+
+ for (typename DoFHandler<dim>::active_cell_iterator
+ cell = scalar_velocity_dof_handler.begin_active();
+ cell != scalar_velocity_dof_handler.end();
+ ++cell)
+ {
+ cell->get_dof_indices (scalar_velocity_cell_dof_indices);
+ scalar_to_vector[cell]->get_dof_indices (cell_dof_indices);
+
+ for (unsigned int i=0; i<scalar_velocity_fe.dofs_per_cell; ++i)
+ scalar_to_vector_velocities[component][scalar_velocity_cell_dof_indices[i]]
+ = cell_dof_indices[fe.component_to_system_index(component, i)];
+ }
+
+ Assert (std::find (scalar_to_vector_velocities[component].begin(),
+ scalar_to_vector_velocities[component].end(),
+ numbers::invalid_unsigned_int)
+ == scalar_to_vector_velocities[component].end(),
+ ExcInternalError());
+ }
+ }
+
// As last action in this function,
// we need to set the vectors
std::cout << " Rebuilding preconditioner..." << std::flush;
- A_preconditioner
+ MatrixTools::create_laplace_matrix (scalar_velocity_dof_handler,
+ QGauss<dim> (degree+1),
+ scalar_velocity_system_matrix);
+ scalar_velocity_hanging_node_constraints.condense (scalar_velocity_system_matrix);
+ {
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (scalar_velocity_dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ boundary_values);
+ Vector<double> dummy (scalar_velocity_dof_handler.n_dofs());
+ MatrixTools::apply_boundary_values (boundary_values,
+ scalar_velocity_system_matrix,
+ dummy,
+ dummy);
+ }
+
+
+ A_preconditioner
= boost::shared_ptr<typename InnerPreconditioner<dim>::type>
(new typename InnerPreconditioner<dim>::type());
A_preconditioner->initialize (system_matrix.block(0,0),