From 41d16b5a8ecf6d6efa0032be1b7de7a23f882d8d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 19 Apr 2011 03:22:02 +0000 Subject: [PATCH] Break out the assembly of the interface term into its own function. This will make it simpler to deal with adaptively refined interfaces. git-svn-id: https://svn.dealii.org/trunk@23610 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-46/step-46.cc | 92 +++++++++++++++++++---------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index 767a479511..7d6639cada 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -73,6 +73,12 @@ class FluidStructureProblem void setup_dofs (); void assemble_system (); + void assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &stokes_fe_face_values, + std::vector > &elasticity_phi, + std::vector > &stokes_phi_grads_u, + std::vector &stokes_phi_p, + FullMatrix &local_interface_matrix) const; void solve (); void output_results (const unsigned int refinement_cycle) const; void refine_mesh (); @@ -94,6 +100,7 @@ class FluidStructureProblem Vector solution; Vector system_rhs; + const double viscosity; const double lambda; const double mu; }; @@ -193,6 +200,7 @@ FluidStructureProblem::FluidStructureProblem (const unsigned int stokes_deg FE_Nothing(), 1, FE_Q(elasticity_degree), dim), dof_handler (triangulation), + viscosity (2), lambda (1), mu (1) { @@ -341,8 +349,6 @@ void FluidStructureProblem::assemble_system () system_matrix=0; system_rhs=0; - const double viscosity = 2; - const QGauss stokes_quadrature(stokes_degree+2); const QGauss elasticity_quadrature(elasticity_degree+2); @@ -493,35 +499,16 @@ void FluidStructureProblem::assemble_system () elasticity_fe_face_values.reinit (cell->neighbor(f), cell->neighbor_of_neighbor(f)); - local_interface_matrix = 0; - for (unsigned int q=0; q normal_vector = stokes_fe_face_values.normal_vector(q); - - for (unsigned int k=0; kneighbor(f)->get_dof_indices (neighbor_dof_indices); - constraints.distribute_local_to_global(local_interface_matrix, - neighbor_dof_indices, - local_dof_indices, - system_matrix); - } + 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, + system_matrix); } else if ((cell->neighbor(f)->level() == cell->level()) && @@ -541,9 +528,50 @@ void FluidStructureProblem::assemble_system () +template +void +FluidStructureProblem::assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &stokes_fe_face_values, + std::vector > &elasticity_phi, + std::vector > &stokes_phi_grads_u, + std::vector &stokes_phi_p, + FullMatrix &local_interface_matrix) const +{ + Assert (stokes_fe_face_values.n_quadrature_points == + elasticity_fe_face_values.n_quadrature_points, + ExcInternalError()); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector displacements (dim+1); + + local_interface_matrix = 0; + for (unsigned int q=0; q normal_vector = stokes_fe_face_values.normal_vector(q); + + for (unsigned int k=0; k -void FluidStructureProblem::solve () +void +FluidStructureProblem::solve () { SparseDirectUMFPACK direct_solver; direct_solver.initialize (system_matrix); -- 2.39.5