From 13cdabe9b4ad9d31cfff164bfdecf3ca709b369c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 19 Apr 2011 21:54:53 +0000 Subject: [PATCH] Implement the missing rest (in particular integrating over refined faces). git-svn-id: https://svn.dealii.org/trunk@23613 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-46/step-46.cc | 85 +++++++++++++++++++++-------- 1 file changed, 62 insertions(+), 23 deletions(-) diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index fa1fba9b13..da3e82cf31 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -377,6 +377,12 @@ void FluidStructureProblem::assemble_system () update_gradients); FEFaceValues elasticity_fe_face_values (elasticity_fe, face_quadrature, update_values); + FESubfaceValues stokes_fe_subface_values (stokes_fe, face_quadrature, + update_JxW_values | + update_normal_vectors | + update_gradients); + FESubfaceValues elasticity_fe_subface_values (elasticity_fe, face_quadrature, + update_values); const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell; @@ -386,7 +392,8 @@ void FluidStructureProblem::assemble_system () stokes_dofs_per_cell); Vector local_rhs; - std::vector local_dof_indices, neighbor_dof_indices; + std::vector local_dof_indices; + std::vector neighbor_dof_indices (stokes_dofs_per_cell); const RightHandSide right_hand_side; @@ -471,8 +478,7 @@ void FluidStructureProblem::assemble_system () } } - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); + local_dof_indices.resize (cell->get_fe().dofs_per_cell); cell->get_dof_indices (local_dof_indices); // local_rhs==0, but need to do @@ -483,50 +489,83 @@ void FluidStructureProblem::assemble_system () system_matrix, system_rhs); // see about face terms - if (cell_is_in_fluid_domain (cell)) - // we are on a fluid cell + if (cell_is_in_solid_domain (cell)) + // we are on a solid cell for (unsigned int f=0; f::faces_per_cell; ++f) if (cell->at_boundary(f) == false) { if ((cell->neighbor(f)->level() == cell->level()) && - (cell->neighbor(f)->has_children() == false)) + (cell->neighbor(f)->has_children() == false) + && + cell_is_in_fluid_domain (cell->neighbor(f))) { - // cells are on same level - if (cell_is_in_fluid_domain (cell->neighbor(f))) - // neighbor is also a fluid cell - continue; - // same size // neighbors; // neighbor is - // solids cell - stokes_fe_face_values.reinit (cell, f); - elasticity_fe_face_values.reinit (cell->neighbor(f), - cell->neighbor_of_neighbor(f)); + // fluid cell + 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_phi_grads_u, stokes_phi_p, local_interface_matrix); - neighbor_dof_indices.resize (elasticity_dofs_per_cell); cell->neighbor(f)->get_dof_indices (neighbor_dof_indices); constraints.distribute_local_to_global(local_interface_matrix, - neighbor_dof_indices, - local_dof_indices, + local_dof_indices, + neighbor_dof_indices, system_matrix); } else if ((cell->neighbor(f)->level() == cell->level()) && (cell->neighbor(f)->has_children() == true)) { - // neighbor has children + // neighbor has children. loop over + // the cells adjacent to the commone + // interface and see which subdomain + // they belong to + for (unsigned int subface=0; subfaceface(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_phi_grads_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); + } } - else + else if (cell->neighbor_is_coarser(f) + && + cell_is_in_fluid_domain(cell->neighbor(f))) { // neighbor is coarser - if (cell->active_fe_index() == cell->neighbor(f)->active_fe_index()) - continue; + 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_phi_grads_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); + } } } @@ -613,7 +652,7 @@ FluidStructureProblem::output_results (const unsigned int refinement_cycle) data_out.add_data_vector (solution, solution_names, DataOut >::type_dof_data, data_component_interpretation); - data_out.build_patches (stokes_fe.degree); + data_out.build_patches (); std::ostringstream filename; filename << "solution-" -- 2.39.5