From 27822def372a875e2e1ca9512b90f83a267c9890 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 23 Feb 2011 15:03:01 +0000 Subject: [PATCH] Make it run with the hp elements. The solution is still wrong, though. git-svn-id: https://svn.dealii.org/trunk@23431 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-46/step-46.cc | 120 ++++++++++++++++++++-------- 1 file changed, 86 insertions(+), 34 deletions(-) diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index f088ac7680..ec51f8edc7 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -86,6 +86,9 @@ class StokesProblem Vector solution; Vector system_rhs; + + const double lambda; + const double mu; }; @@ -182,7 +185,9 @@ StokesProblem::StokesProblem (const unsigned int stokes_degree, elasticity_fe (FE_Nothing(), dim, FE_Nothing(), 1, FE_Q(elasticity_degree), dim), - dof_handler (triangulation) + dof_handler (triangulation), + lambda (1), + mu (1) { fe_collection.push_back (stokes_fe); fe_collection.push_back (elasticity_fe); @@ -207,9 +212,8 @@ void StokesProblem::setup_dofs () && (cell->center()[1] > -0.5))) cell->set_active_fe_index (0); -//TODO: doesn't work right now due to domination issues - // else - // cell->set_active_fe_index (1); + else + cell->set_active_fe_index (1); dof_handler.distribute_dofs (fe_collection); DoFRenumbering::Cuthill_McKee (dof_handler); @@ -259,10 +263,10 @@ void StokesProblem::assemble_system () system_matrix=0; system_rhs=0; - QGauss stokes_quadrature(stokes_degree+2); - QGauss elasticity_quadrature(elasticity_degree+2); + const QGauss stokes_quadrature(stokes_degree+2); + const QGauss elasticity_quadrature(elasticity_degree+2); - hp::QCollection q_collection; + hp::QCollection q_collection; q_collection.push_back (stokes_quadrature); q_collection.push_back (elasticity_quadrature); @@ -272,23 +276,26 @@ void StokesProblem::assemble_system () update_JxW_values | update_gradients); - const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; - - const unsigned int n_q_points = stokes_quadrature.size(); + const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell; + const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell; - FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); - Vector local_rhs (dofs_per_cell); + FullMatrix local_matrix; + Vector local_rhs; - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices; const RightHandSide right_hand_side; const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector displacements (dim+1); + + std::vector > stokes_phi_grads_u (stokes_dofs_per_cell); + std::vector stokes_div_phi_u (stokes_dofs_per_cell); + std::vector stokes_phi_p (stokes_dofs_per_cell); - std::vector > phi_grads_u (dofs_per_cell); - std::vector div_phi_u (dofs_per_cell); - std::vector phi_p (dofs_per_cell); + std::vector > elasticity_phi_grad (elasticity_dofs_per_cell); + std::vector elasticity_phi_div (elasticity_dofs_per_cell); typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -299,31 +306,69 @@ void StokesProblem::assemble_system () const FEValues &fe_values = hp_fe_values.get_present_fe_values(); - local_matrix = 0; - local_rhs = 0; + local_matrix.reinit (cell->get_fe().dofs_per_cell, + cell->get_fe().dofs_per_cell); + local_rhs.reinit (cell->get_fe().dofs_per_cell); - for (unsigned int q=0; qactive_fe_index() == 0) { - for (unsigned int k=0; kget_fe().dofs_per_cell; + Assert (dofs_per_cell == stokes_dofs_per_cell, + ExcInternalError()); + + for (unsigned int q=0; qget_fe().dofs_per_cell; + Assert (dofs_per_cell == elasticity_dofs_per_cell, + ExcInternalError()); + for (unsigned int q=0; qget_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); // local_rhs==0, but need to do @@ -394,10 +439,17 @@ StokesProblem::refine_mesh () { Vector estimated_error_per_cell (triangulation.n_active_cells()); + const QGauss stokes_quadrature(stokes_degree+2); + const QGauss elasticity_quadrature(elasticity_degree+2); + + hp::QCollection face_q_collection; + face_q_collection.push_back (stokes_quadrature); + face_q_collection.push_back (elasticity_quadrature); + std::vector component_mask (dim+1+dim, false); component_mask[dim] = true; KellyErrorEstimator::estimate (dof_handler, - QGauss(stokes_degree+1), + face_q_collection, typename FunctionMap::type(), solution, estimated_error_per_cell, -- 2.39.5