From: Wolfgang Bangerth Date: Fri, 4 Dec 2020 03:25:58 +0000 (-0700) Subject: Fix step-46. X-Git-Tag: v9.3.0-rc1~803^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=51d97d5ae791293d8c2e9d8f1519b2b754169d38;p=dealii.git Fix step-46. --- diff --git a/examples/step-46/step-46.cc b/examples/step-46/step-46.cc index f33a9890c4..ac68526915 100644 --- a/examples/step-46/step-46.cc +++ b/examples/step-46/step-46.cc @@ -367,35 +367,37 @@ namespace Step46 for (const auto &cell : dof_handler.active_cell_iterators()) if (cell_is_in_fluid_domain(cell)) for (const auto face_no : cell->face_indices()) - { - bool face_is_on_interface = false; - - if ((cell->neighbor(face_no)->has_children() == false) && - (cell_is_in_solid_domain(cell->neighbor(face_no)))) - face_is_on_interface = true; - else if (cell->neighbor(face_no)->has_children() == true) - { - for (unsigned int sf = 0; - sf < cell->face(face_no)->n_children(); - ++sf) - if (cell_is_in_solid_domain( - cell->neighbor_child_on_subface(face_no, sf))) - { - face_is_on_interface = true; - break; - } - } - - if (face_is_on_interface) - { - cell->face(face_no)->get_dof_indices(local_face_dof_indices, - 0); - for (unsigned int i = 0; i < local_face_dof_indices.size(); - ++i) - if (stokes_fe.face_system_to_component_index(i).first < dim) - constraints.add_line(local_face_dof_indices[i]); - } - } + if (cell->face(face_no)->at_boundary() == false) + { + bool face_is_on_interface = false; + + if ((cell->neighbor(face_no)->has_children() == false) && + (cell_is_in_solid_domain(cell->neighbor(face_no)))) + face_is_on_interface = true; + else if (cell->neighbor(face_no)->has_children() == true) + { + for (unsigned int sf = 0; + sf < cell->face(face_no)->n_children(); + ++sf) + if (cell_is_in_solid_domain( + cell->neighbor_child_on_subface(face_no, sf))) + { + face_is_on_interface = true; + break; + } + } + + if (face_is_on_interface) + { + cell->face(face_no)->get_dof_indices(local_face_dof_indices, + 0); + for (unsigned int i = 0; i < local_face_dof_indices.size(); + ++i) + if (stokes_fe.face_system_to_component_index(i).first < + dim) + constraints.add_line(local_face_dof_indices[i]); + } + } } // At the end of all this, we can declare to the constraints object that @@ -641,118 +643,121 @@ namespace Step46 // domain. Let's start with these conditions: if (cell_is_in_solid_domain(cell)) for (const auto f : cell->face_indices()) - { - // At this point we know that the current cell is a candidate - // for integration and that a neighbor behind face - // f exists. There are now three possibilities: - // - // - The neighbor is at the same refinement level and has no - // children. - // - The neighbor has children. - // - The neighbor is coarser. - // - // In all three cases, we are only interested in it if it is - // part of the fluid subdomain. So let us start with the first - // and simplest case: if the neighbor is at the same level, - // has no children, and is a fluid cell, then the two cells - // share a boundary that is part of the interface along which - // we want to integrate interface terms. All we have to do is - // initialize two FEFaceValues object with the current face - // and the face of the neighboring cell (note how we find out - // which face of the neighboring cell borders on the current - // cell) and pass things off to the function that evaluates - // the interface terms (the third through fifth arguments to - // this function provide it with scratch arrays). The result - // is then again copied into the global matrix, using a - // function that knows that the DoF indices of rows and - // columns of the local matrix result from different cells: - if ((cell->neighbor(f)->level() == cell->level()) && - (cell->neighbor(f)->has_children() == false) && - cell_is_in_fluid_domain(cell->neighbor(f))) - { - elasticity_fe_face_values.reinit(cell, f); - stokes_fe_face_values.reinit(cell->neighbor(f), - cell->neighbor_of_neighbor(f)); - - assemble_interface_term(elasticity_fe_face_values, - stokes_fe_face_values, - elasticity_phi, - stokes_symgrad_phi_u, - stokes_phi_p, - local_interface_matrix); - - cell->neighbor(f)->get_dof_indices(neighbor_dof_indices); - constraints.distribute_local_to_global(local_interface_matrix, - local_dof_indices, - neighbor_dof_indices, - system_matrix); - } - - // The second case is if the neighbor has further children. In - // that case, we have to loop over all the children of the - // neighbor to see if they are part of the fluid subdomain. If - // they are, then we integrate over the common interface, - // which is a face for the neighbor and a subface of the - // current cell, requiring us to use an FEFaceValues for the - // neighbor and an FESubfaceValues for the current cell: - else if ((cell->neighbor(f)->level() == cell->level()) && - (cell->neighbor(f)->has_children() == true)) - { - for (unsigned int subface = 0; - subface < cell->face(f)->n_children(); - ++subface) - if (cell_is_in_fluid_domain( - cell->neighbor_child_on_subface(f, subface))) - { - elasticity_fe_subface_values.reinit(cell, f, subface); - stokes_fe_face_values.reinit( - cell->neighbor_child_on_subface(f, subface), - cell->neighbor_of_neighbor(f)); - - assemble_interface_term(elasticity_fe_subface_values, - stokes_fe_face_values, - elasticity_phi, - stokes_symgrad_phi_u, - stokes_phi_p, - local_interface_matrix); - - cell->neighbor_child_on_subface(f, subface) - ->get_dof_indices(neighbor_dof_indices); - constraints.distribute_local_to_global( - local_interface_matrix, - local_dof_indices, - neighbor_dof_indices, - system_matrix); - } - } - - // The last option is that the neighbor is coarser. In that - // case we have to use an FESubfaceValues object for the - // neighbor and a FEFaceValues for the current cell; the rest - // is the same as before: - else if (cell->neighbor_is_coarser(f) && - cell_is_in_fluid_domain(cell->neighbor(f))) - { - elasticity_fe_face_values.reinit(cell, f); - stokes_fe_subface_values.reinit( - cell->neighbor(f), - cell->neighbor_of_coarser_neighbor(f).first, - cell->neighbor_of_coarser_neighbor(f).second); - - assemble_interface_term(elasticity_fe_face_values, - stokes_fe_subface_values, - elasticity_phi, - stokes_symgrad_phi_u, - stokes_phi_p, - local_interface_matrix); - - cell->neighbor(f)->get_dof_indices(neighbor_dof_indices); - constraints.distribute_local_to_global(local_interface_matrix, - local_dof_indices, - neighbor_dof_indices, - system_matrix); - } - } + if (cell->face(f)->at_boundary() == false) + { + // At this point we know that the current cell is a candidate + // for integration and that a neighbor behind face + // f exists. There are now three possibilities: + // + // - The neighbor is at the same refinement level and has no + // children. + // - The neighbor has children. + // - The neighbor is coarser. + // + // In all three cases, we are only interested in it if it is + // part of the fluid subdomain. So let us start with the first + // and simplest case: if the neighbor is at the same level, + // has no children, and is a fluid cell, then the two cells + // share a boundary that is part of the interface along which + // we want to integrate interface terms. All we have to do is + // initialize two FEFaceValues object with the current face + // and the face of the neighboring cell (note how we find out + // which face of the neighboring cell borders on the current + // cell) and pass things off to the function that evaluates + // the interface terms (the third through fifth arguments to + // this function provide it with scratch arrays). The result + // is then again copied into the global matrix, using a + // function that knows that the DoF indices of rows and + // columns of the local matrix result from different cells: + if ((cell->neighbor(f)->level() == cell->level()) && + (cell->neighbor(f)->has_children() == false) && + cell_is_in_fluid_domain(cell->neighbor(f))) + { + elasticity_fe_face_values.reinit(cell, f); + stokes_fe_face_values.reinit(cell->neighbor(f), + cell->neighbor_of_neighbor(f)); + + assemble_interface_term(elasticity_fe_face_values, + stokes_fe_face_values, + elasticity_phi, + stokes_symgrad_phi_u, + stokes_phi_p, + local_interface_matrix); + + cell->neighbor(f)->get_dof_indices(neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + + // The second case is if the neighbor has further children. In + // that case, we have to loop over all the children of the + // neighbor to see if they are part of the fluid subdomain. If + // they are, then we integrate over the common interface, + // which is a face for the neighbor and a subface of the + // current cell, requiring us to use an FEFaceValues for the + // neighbor and an FESubfaceValues for the current cell: + else if ((cell->neighbor(f)->level() == cell->level()) && + (cell->neighbor(f)->has_children() == true)) + { + for (unsigned int subface = 0; + subface < cell->face(f)->n_children(); + ++subface) + if (cell_is_in_fluid_domain( + cell->neighbor_child_on_subface(f, subface))) + { + elasticity_fe_subface_values.reinit(cell, f, subface); + stokes_fe_face_values.reinit( + cell->neighbor_child_on_subface(f, subface), + cell->neighbor_of_neighbor(f)); + + assemble_interface_term(elasticity_fe_subface_values, + stokes_fe_face_values, + elasticity_phi, + stokes_symgrad_phi_u, + stokes_phi_p, + local_interface_matrix); + + cell->neighbor_child_on_subface(f, subface) + ->get_dof_indices(neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + } + + // The last option is that the neighbor is coarser. In that + // case we have to use an FESubfaceValues object for the + // neighbor and a FEFaceValues for the current cell; the rest + // is the same as before: + else if (cell->neighbor_is_coarser(f) && + cell_is_in_fluid_domain(cell->neighbor(f))) + { + elasticity_fe_face_values.reinit(cell, f); + stokes_fe_subface_values.reinit( + cell->neighbor(f), + cell->neighbor_of_coarser_neighbor(f).first, + cell->neighbor_of_coarser_neighbor(f).second); + + assemble_interface_term(elasticity_fe_face_values, + stokes_fe_subface_values, + elasticity_phi, + stokes_symgrad_phi_u, + stokes_phi_p, + local_interface_matrix); + + cell->neighbor(f)->get_dof_indices(neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + } } }