From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 19 Apr 2011 14:26:54 +0000 (+0000)
Subject: Split out setting the subdomain indicators. Use both velocity and displacement to... 
X-Git-Tag: v8.0.0~4145
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8e7fe293b44d44882d9d8b7ef2547596f8040dd7;p=dealii.git

Split out setting the subdomain indicators. Use both velocity and displacement to refine the mesh.

git-svn-id: https://svn.dealii.org/trunk@23612 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc
index 7d6639cada..fa1fba9b13 100644
--- a/deal.II/examples/step-46/step-46.cc
+++ b/deal.II/examples/step-46/step-46.cc
@@ -71,6 +71,7 @@ class FluidStructureProblem
     cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell);
 
 
+    void setup_subdomains ();
     void setup_dofs ();
     void assemble_system ();
     void assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
@@ -235,20 +236,6 @@ void FluidStructureProblem<dim>::setup_dofs ()
 {
   system_matrix.clear ();
 
-  for (typename hp::DoFHandler<dim>::active_cell_iterator
-	 cell = dof_handler.begin_active();
-       cell != dof_handler.end(); ++cell)
-    if (((std::fabs(cell->center()[0]) < 0.25)
-	 &&
-	 (cell->center()[1] > 0.5))
-	||
-	((std::fabs(cell->center()[0]) >= 0.25)
-	 &&
-	 (cell->center()[1] > -0.5)))
-      cell->set_active_fe_index (0);
-    else
-      cell->set_active_fe_index (1);
-
   dof_handler.distribute_dofs (fe_collection);
   DoFRenumbering::Cuthill_McKee (dof_handler);
 
@@ -342,6 +329,25 @@ void FluidStructureProblem<dim>::setup_dofs ()
 }
 
 
+template <int dim>
+void
+FluidStructureProblem<dim>::setup_subdomains ()
+{
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    if (((std::fabs(cell->center()[0]) < 0.25)
+         &&
+         (cell->center()[1] > 0.5))
+	||
+	((std::fabs(cell->center()[0]) >= 0.25)
+	 &&
+	 (cell->center()[1] > -0.5)))
+      cell->set_active_fe_index (0);
+    else
+      cell->set_active_fe_index (1);
+}
+
 
 template <int dim>
 void FluidStructureProblem<dim>::assemble_system ()
@@ -624,6 +630,8 @@ template <int dim>
 void
 FluidStructureProblem<dim>::refine_mesh ()
 {
+  Vector<float> stokes_estimated_error_per_cell (triangulation.n_active_cells());
+  Vector<float> elasticity_estimated_error_per_cell (triangulation.n_active_cells());
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
   const QGauss<dim-1> stokes_quadrature(stokes_degree+2);
@@ -633,15 +641,30 @@ FluidStructureProblem<dim>::refine_mesh ()
   face_q_collection.push_back (stokes_quadrature);
   face_q_collection.push_back (elasticity_quadrature);
 
-  std::vector<bool> component_mask (dim+1+dim, false);
-  component_mask[dim] = true;
+  std::vector<bool> stokes_component_mask (dim+1+dim, false);
+  for (unsigned int d=0; d<dim; ++d)
+    stokes_component_mask[d] = true;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                      face_q_collection,
+                                      typename FunctionMap<dim>::type(),
+                                      solution,
+                                      stokes_estimated_error_per_cell,
+                                      stokes_component_mask);
+
+  std::vector<bool> elasticity_component_mask (dim+1+dim, false);
+  for (unsigned int d=0; d<dim; ++d)
+    elasticity_component_mask[dim+1+d] = true;
   KellyErrorEstimator<dim>::estimate (dof_handler,
                                       face_q_collection,
                                       typename FunctionMap<dim>::type(),
                                       solution,
-                                      estimated_error_per_cell,
-                                      component_mask);
+                                      elasticity_estimated_error_per_cell,
+                                      elasticity_component_mask);
 
+  stokes_estimated_error_per_cell /= stokes_estimated_error_per_cell.linfty_norm();
+  elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.linfty_norm();
+  estimated_error_per_cell += stokes_estimated_error_per_cell;
+  estimated_error_per_cell += elasticity_estimated_error_per_cell;
   GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                    estimated_error_per_cell,
                                                    0.3, 0.0);
@@ -672,6 +695,7 @@ void FluidStructureProblem<dim>::run ()
       if (refinement_cycle > 0)
         refine_mesh ();
 
+      setup_subdomains ();
       setup_dofs ();
 
       std::cout << "   Assembling..." << std::endl << std::flush;