From a594840381186118a96574cb67af4bb7bcdb7540 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 5 Aug 2011 22:01:26 +0000 Subject: [PATCH] Use a different formulation where we use a different enriched space that is only defined on cells that are intersected and where the additional basis functions are zero at the vertices of cells that are intersected. git-svn-id: https://svn.dealii.org/trunk@24026 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-47/step-47.cc | 122 ++++++++-------------------- 1 file changed, 34 insertions(+), 88 deletions(-) diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index df97516c7a..5f5b3a8bd2 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -229,48 +229,13 @@ interface_intersects_cell (const typename Triangulation::cell_iterator &cel template void LaplaceProblem::setup_system () { - // 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); - } + if (interface_intersects_cell(cell) == false) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); dof_handler.distribute_dofs (fe_collection); @@ -283,26 +248,6 @@ void LaplaceProblem::setup_system () // DoFTools::make_hanging_node_constraints (dof_handler, // constraints); - // 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) - constraints.add_line (cell->vertex_dof_index(v,1,cell->active_fe_index())); - //TODO: component 1 must satisfy zero boundary conditions constraints.close(); @@ -317,6 +262,7 @@ void LaplaceProblem::setup_system () system_matrix.reinit (sparsity_pattern); } + template void LaplaceProblem::assemble_system () { @@ -377,32 +323,21 @@ void LaplaceProblem::assemble_system () } else { +//TODO: verify that the order of support points equals the order of vertices of the cells, as we use below +//TODO: remove update_support_points and friends, since they aren't implemented anyway Assert (cell->active_fe_index() == 1, ExcInternalError()); + Assert (interface_intersects_cell(cell) == true, ExcInternalError()); - std::auto_ptr > custom_fe_values; + std::vector level_set_values (GeometryInfo::vertices_per_cell); + for (unsigned int v=0; v::vertices_per_cell; ++v) + level_set_values[v] = level_set (cell->vertex(v)); - if (interface_intersects_cell(cell) == false) - enriched_fe_values.reinit (cell); - else - { - std::vector level_set_values (GeometryInfo::vertices_per_cell); - for (unsigned int v=0; v::vertices_per_cell; ++v) - level_set_values[v] = level_set (cell->vertex(v)); - - custom_fe_values - .reset(new FEValues (fe_collection[1], - compute_quadrature(quadrature_formula, cell, - level_set_values).second, - update_values | update_gradients | - update_quadrature_points | update_JxW_values)); - custom_fe_values->reinit (cell); - } - - FEValues &this_fe_values = ((interface_intersects_cell(cell) == false) - ? - enriched_fe_values - : - *custom_fe_values); + FEValues this_fe_values (fe_collection[1], + compute_quadrature(quadrature_formula, cell, + level_set_values).second, + update_values | update_gradients | + update_quadrature_points | update_JxW_values | + update_support_points); coefficient_values.resize (this_fe_values.n_quadrature_points); coefficient.value_list (this_fe_values.get_quadrature_points(), @@ -420,8 +355,11 @@ void LaplaceProblem::assemble_system () this_fe_values.JxW(q_point)); else cell_matrix(i,j) += (coefficient_values[q_point] * - this_fe_values.shape_grad(i,q_point) * - (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + this_fe_values.shape_grad(i,q_point) + * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* this_fe_values.shape_grad(j,q_point) + grad_level_set(this_fe_values.quadrature_point(q_point)) * @@ -438,7 +376,9 @@ void LaplaceProblem::assemble_system () for (unsigned int j=0; jget_fe().system_to_component_index(j).first == 0) cell_matrix(i,j) += (coefficient_values[q_point] * - (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* this_fe_values.shape_grad(i,q_point) + grad_level_set(this_fe_values.quadrature_point(q_point)) * @@ -448,13 +388,17 @@ void LaplaceProblem::assemble_system () this_fe_values.JxW(q_point)); else cell_matrix(i,j) += (coefficient_values[q_point] * - (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* this_fe_values.shape_grad(i,q_point) + grad_level_set(this_fe_values.quadrature_point(q_point)) * sign(level_set(this_fe_values.quadrature_point(q_point))) * this_fe_values.shape_value(i,q_point)) * - (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(j).second))))* this_fe_values.shape_grad(j,q_point) + grad_level_set(this_fe_values.quadrature_point(q_point)) * @@ -462,7 +406,9 @@ void LaplaceProblem::assemble_system () this_fe_values.shape_value(j,q_point)) * this_fe_values.JxW(q_point)); - cell_rhs(i) += (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) * + cell_rhs(i) += ((std::fabs(level_set(this_fe_values.quadrature_point(q_point))) + - + std::fabs(level_set(cell->vertex(cell->get_fe().system_to_component_index(i).second))))* this_fe_values.shape_value(i,q_point) * 1.0 * this_fe_values.JxW(q_point)); -- 2.39.5