From: bangerth Date: Tue, 2 Aug 2011 20:11:42 +0000 (+0000) Subject: Implement the first few steps of what's necessary here. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=96840a55c9e2764f34df80e60fbbe7170135c3b5;p=dealii-svn.git Implement the first few steps of what's necessary here. git-svn-id: https://svn.dealii.org/trunk@23997 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index f3165c878e..d3b359511d 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -28,6 +28,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -60,24 +63,26 @@ class LaplaceProblem void run (); private: + bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; + void setup_system (); void assemble_system (); void solve (); void refine_grid (); void output_results (const unsigned int cycle) const; - Triangulation triangulation; + Triangulation triangulation; - DoFHandler dof_handler; - FE_Q fe; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; - ConstraintMatrix hanging_node_constraints; + ConstraintMatrix hanging_node_constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; + Vector solution; + Vector system_rhs; }; @@ -139,9 +144,13 @@ void Coefficient::value_list (const std::vector > &points, template LaplaceProblem::LaplaceProblem () : - dof_handler (triangulation), - fe (2) -{} + dof_handler (triangulation) +{ + fe_collection.push_back (FESystem (FE_Q(1), 1, + FE_Nothing(), 1)); + fe_collection.push_back (FESystem (FE_Q(1), 1, + FE_Q(1), 1)); +} @@ -153,10 +162,63 @@ LaplaceProblem::~LaplaceProblem () +template +bool +LaplaceProblem:: +interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const +{ + return false; +} + + + template void LaplaceProblem::setup_system () { - dof_handler.distribute_dofs (fe); + // decide which element to use + // where. to do so, we need to know + // which elements are intersected + // by the interface, or are at + // least adjacent. to this end, in + // a first step, loop over all + // cells and record which vertices + // are on cells that are + // intersected + std::vector vertex_is_on_intersected_cell (triangulation.n_vertices(), + false); + for (typename Triangulation::cell_iterator cell + = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (interface_intersects_cell(cell)) + for (unsigned int v=0; v::vertices_per_cell; ++v) + vertex_is_on_intersected_cell[cell->vertex_index(v)] = true; + + // now loop over all cells + // again. if one of the vertices of + // a given cell is part of a cell + // that is intersected, then we + // need to use the enriched space + // there. otherwise, use the normal + // space + for (typename hp::DoFHandler::cell_iterator cell + = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + { + bool use_enriched_space = false; + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertex_is_on_intersected_cell[cell->vertex_index(v)] == true) + { + use_enriched_space = true; + break; + } + + if (use_enriched_space == false) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + + dof_handler.distribute_dofs (fe_collection); solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); @@ -165,9 +227,30 @@ void LaplaceProblem::setup_system () hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); - hanging_node_constraints.close (); + // now constrain those enriched + // DoFs that are on cells that are + // not intersected but that are + // adjacent to cells that are + for (typename hp::DoFHandler::cell_iterator cell + = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if ((cell->active_fe_index() == 1) + && + (interface_intersects_cell(cell) == false)) + // we are on an enriched cell + // but it isn't + // intersected. see which + // vertices are not part of + // intersected cells and + // constrain these DoFs + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertex_is_on_intersected_cell[cell->vertex_index(v)] == false) + hanging_node_constraints.add_line (cell->vertex_dof_index(v,1)); + hanging_node_constraints.close(); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); @@ -184,32 +267,35 @@ void LaplaceProblem::assemble_system () { const QGauss quadrature_formula(3); - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); + FEValues plain_fe_values (fe_collection[0], 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(); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + FullMatrix cell_matrix; + Vector cell_rhs; - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices; const Coefficient coefficient; std::vector coefficient_values (n_q_points); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + cell_rhs.reinit (dofs_per_cell); + cell_matrix = 0; cell_rhs = 0; - fe_values.reinit (cell); + plain_fe_values.reinit (cell); - coefficient.value_list (fe_values.get_quadrature_points(), + coefficient.value_list (plain_fe_values.get_quadrature_points(), coefficient_values); for (unsigned int q_point=0; q_point::assemble_system () { for (unsigned int j=0; jget_dof_indices (local_dof_indices); - for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, - ZeroFunction(), + ZeroFunction(2), boundary_values); MatrixTools::apply_boundary_values (boundary_values, system_matrix, @@ -349,7 +427,7 @@ void LaplaceProblem::run () DataOutBase::EpsFlags eps_flags; eps_flags.z_scaling = 4; - DataOut data_out; + DataOut > data_out; data_out.set_flags (eps_flags); data_out.attach_dof_handler (dof_handler);