[cig-commits] commit: Change it from lid-driven to pure shear.

[Mercurial]

changeset:   1:99033ff21ab4
user:        Walter Landry <wlandry at caltech.edu>
date:        Sat Feb 19 11:35:27 2011 -0800
files:       madds-1.cc
description:
Change it from lid-driven to pure shear.


diff -r 62413b0e7b38 -r 99033ff21ab4 madds-1.cc
--- a/madds-1.cc	Fri Feb 18 16:50:04 2011 -0800
+++ b/madds-1.cc	Sat Feb 19 11:35:27 2011 -0800
@@ -305,13 +305,11 @@ TopBoundaryValues<dim>::value (const Poi
 {
   Assert (component < this->n_components,
           ExcIndexRange (component, 0, this->n_components));
-
-  if (component == 0)
+  
+  if (component == 1 && p[1]>0.999999)
     return 1;
   return 0;
 }
-
-
 
 template <int dim>
 void
@@ -320,6 +318,44 @@ TopBoundaryValues<dim>::vector_value (co
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = TopBoundaryValues<dim>::value (p, c);
+}
+
+template <int dim>
+class RightBoundaryValues : public Function<dim>
+{
+  public:
+    RightBoundaryValues ()
+		    :
+		    Function<dim>(dim+1)
+      {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &value) const;
+};
+
+template <int dim>
+double
+RightBoundaryValues<dim>::value (const Point<dim>  &p,
+			       const unsigned int component) const
+{
+  Assert (component < this->n_components,
+          ExcIndexRange (component, 0, this->n_components));
+
+  if (component == 0 && p[0]>0.999999)
+    return -1;
+  return 0;
+}
+
+template <int dim>
+void
+RightBoundaryValues<dim>::vector_value (const Point<dim> &p,
+				      Vector<double>   &values) const
+{
+  for (unsigned int c=0; c<this->n_components; ++c)
+    values(c) = RightBoundaryValues<dim>::value (p, c);
 }
 
 
@@ -444,16 +480,27 @@ void StokesProblem<dim>::setup_dofs ()
     constraints.clear ();
     std::vector<bool> component_mask (dim+1, true);
     component_mask[dim] = false;
+
+    component_mask[dim-2] = false;
     VectorTools::interpolate_boundary_values (dof_handler,
-					      1,
-					      TopBoundaryValues<dim>(),
-					      constraints,
-					      component_mask);
+        				      1,
+        				      TopBoundaryValues<dim>(),
+        				      constraints,
+        				      component_mask);
+    component_mask[dim-1] = false;
+    component_mask[dim-2] = true;
     VectorTools::interpolate_boundary_values (dof_handler,
-					      0,
-					      ZeroFunction<dim>(dim+1),
-					      constraints,
-					      component_mask);
+        				      2,
+        				      RightBoundaryValues<dim>(),
+        				      constraints,
+        				      component_mask);
+
+    // component_mask[dim-1] = true;
+    // VectorTools::interpolate_boundary_values (dof_handler,
+    //     				      0,
+    //     				      ZeroFunction<dim>(dim+1),
+    //     				      constraints,
+    //     				      component_mask);
     DoFTools::make_hanging_node_constraints (dof_handler,
 					     constraints);
   }
@@ -743,8 +790,17 @@ void StokesProblem<dim>::run ()
 	 cell = triangulation.begin_active();
        cell != triangulation.end(); ++cell)
     for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-      if (triangulation.begin_active()->face(f)->center()[dim-1] == 1)
-	cell->face(f)->set_all_boundary_indicators(1);
+      if (triangulation.begin_active()->face(f)->center()[dim-1] == 1
+          || triangulation.begin_active()->face(f)->center()[dim-1] == 0)
+        {
+          cell->face(f)->set_all_boundary_indicators(1);
+        }
+      else if (triangulation.begin_active()->face(f)->center()[dim-2] == 1
+               || triangulation.begin_active()->face(f)->center()[dim-2] == 0)
+        {
+          cell->face(f)->set_all_boundary_indicators(2);
+        }
+
   triangulation.refine_global (parameters.initial_refinement);
 
   for (unsigned int refinement_cycle = 0;