Commits

Garth Wells committed 55bf4b0 Merge

Merge remote-tracking branch 'origin/garth/fixes-snes-vi'

  • Participants
  • Parent commits 83d082c, 33094b6

Comments (0)

Files changed (7)

cmake/scripts/download-demo-data

 }
 
 download_mesh demo/documented/eigenvalue                     box_with_dent.xml.gz
-download_mesh demo/undocumented/meshfunction                 unitsquare_2_2.xml.gz
-download_mesh demo/undocumented/meshfunction                 unitsquare_2_2_subdomains.xml.gz
 download_mesh demo/undocumented/mixed-poisson-sphere         sphere_16.xml.gz
 download_mesh demo/undocumented/advection-diffusion          dolfin_fine.xml.gz
 download_mesh demo/undocumented/advection-diffusion          dolfin_fine_velocity.xml.gz
 download_mesh demo/undocumented/advection-diffusion          dolfin_fine_subdomains.xml.gz
+download_mesh demo/undocumented/contact-vi-snes              circle_yplane.xml.gz
+download_mesh demo/undocumented/contact-vi-tao               circle_yplane.xml.gz
 download_mesh demo/undocumented/multistage-solver            dolfin_fine.xml.gz
 download_mesh demo/undocumented/multistage-solver            dolfin_fine_velocity.xml.gz
 download_mesh demo/undocumented/multistage-solver            dolfin_fine_subdomains.xml.gz
 download_mesh demo/documented/tensor-weighted-poisson        unitsquare_32_32_c11.xml.gz
 download_mesh demo/documented/stokes-taylor-hood/            dolfin_fine.xml.gz
 download_mesh demo/documented/stokes-taylor-hood             dolfin_fine_subdomains.xml.gz
-download_mesh demo/undocumented/stokes-taylor-hood           dolfin_fine.xml.gz
 download_mesh demo/undocumented/elastodynamics               dolfin_fine.xml.gz
 download_mesh demo/undocumented/auto-adaptive-navier-stokes  channel_with_flap.xml.gz
 download_mesh demo/documented/stokes-mini                    dolfin_fine.xml.gz

demo/undocumented/contact-vi-snes/circle_yplane.geo

+// Gmsh input file for unit radius circle with mesh symmetry line on y-axis
+cl__1 = 0.05;
+Point(1) = {0, 1, 0, cl__1};
+Point(2) = {1, 0, 0, cl__1};
+Point(3) = {0, -1, 0, cl__1};
+Point(4) = {-1, 0, 0, cl__1};
+Point(6) = {0, 0, 0, cl__1};
+Circle(1) = {1, 6, 4};
+Circle(2) = {4, 6, 3};
+Circle(3) = {3, 6, 2};
+Circle(4) = {2, 6, 1};
+Line Loop(7) = {1, 2, 3, 4};
+Plane Surface(7) = {7};
+Line(8) = {1, 3};
+Line{8} In Surface {7};

demo/undocumented/contact-vi-snes/cpp/main.cpp

 
 using namespace dolfin;
 
+// Sub domain for symmetry condition
+class SymmetryLine : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  { return (std::abs(x[0]) < DOLFIN_EPS); }
+};
+
+// Lower bound for displacement
+class LowerBound : public Expression
+{
+public:
+  LowerBound() : Expression(2) {}
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    const double xmin = -1.0 - DOLFIN_EPS;
+    const double ymin = -1.0;
+    values[0] = xmin - x[0];
+    values[1] = ymin - x[1];
+  }
+};
+
+// Upper bound for displacement
+class UpperBound : public Expression
+{
+public:
+  UpperBound() : Expression(2) {}
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    const double xmax = 1.0 + DOLFIN_EPS;
+    const double ymax = 2.0;
+    values[0] = xmax - x[0];
+    values[1] = ymax - x[1];
+  }
+};
+
 int main()
 {
 #ifdef HAS_PETSC
-#ifdef HAS_CGAL
-    // Sub domain for symmetry condition
-    class SymmetryLine : public SubDomain
-    {
-        bool inside(const Array<double>& x, bool on_boundary) const
-        {
-        return (std::abs(x[0]) < DOLFIN_EPS);
-        }
-    };
-    // Lower bound for displacement
-    class LowerBound : public Expression
-    {
-    public:
-
-      LowerBound() : Expression(2) {}
-
-      void eval(Array<double>& values, const Array<double>& x) const
-      {
-        double xmin = -1.-DOLFIN_EPS;
-        double ymin = -1.;
-        values[0] = xmin-x[0];
-        values[1] = ymin-x[1];
-      }
-
-    };
-
-    // Upper bound for displacement
-    class UpperBound : public Expression
-    {
-    public:
-
-      UpperBound() : Expression(2) {}
-
-      void eval(Array<double>& values, const Array<double>& x) const
-      {
-        double xmax = 1.+DOLFIN_EPS;
-        double ymax = 2.;
-        values[0] = xmax-x[0];
-        values[1] = ymax-x[1];
-      }
-
-    };
 
   // Read mesh and create function space
-#ifdef HAS_CGAL
-  Circle circle(0, 0, 1);
-  Mesh   mesh(circle,30);
-#else
-  UnitCircleMesh mesh(30);
-#endif
+  Mesh mesh("../circle_yplane.xml.gz");
+
+  // Create function space
   HyperElasticity::FunctionSpace V(mesh);
 
   // Create Dirichlet boundary conditions
   const double E  = 10.0;
   const double nu = 0.3;
   Constant mu(E/(2*(1 + nu)));
-  Constant lambda(E*nu/((1 + nu)*(1 - 2*nu)));
+  Constant lambda(E*nu/((1.0 + nu)*(1.0 - 2.0*nu)));
 
   // Create (linear) form defining (nonlinear) variational problem
   HyperElasticity::ResidualForm F(V);
   HyperElasticity::JacobianForm J(V, V);
   J.mu = mu; J.lmbda = lambda; J.u = u;
 
-  // Interpolate expression for Upper bound
+  // Interpolate expression for upper bound
   UpperBound umax_exp;
   Function umax(V);
   umax.interpolate(umax_exp);
 
-  // Interpolate expression for Lower bound
+  // Interpolate expression for lower bound
   LowerBound umin_exp;
   Function umin(V);
   umin.interpolate(umin_exp);
 
   // Set up the non-linear solver
   NonlinearVariationalSolver solver(problem);
-  solver.parameters["nonlinear_solver"]="snes";
-  solver.parameters["linear_solver"]="lu";
-  solver.parameters("snes_solver")["maximum_iterations"]=20;
-  solver.parameters("snes_solver")["report"]=true;
-  solver.parameters("snes_solver")["error_on_nonconvergence"]=false;
-  //info(solver.parameters,true);
+  solver.parameters["nonlinear_solver"] = "snes";
+  solver.parameters["linear_solver"] = "lu";
+  solver.parameters("snes_solver")["maximum_iterations"] = 20;
+  solver.parameters("snes_solver")["report"] = true;
+  solver.parameters("snes_solver")["error_on_nonconvergence"] = false;
 
   // Solve the problems
   std::pair<std::size_t, bool> out;
   out = solver.solve(umin,umax);
 
-  // Check for convergence. Convergence is one modifies the loading and the mesh size.
+  // Check for convergence. Convergence is one modifies the loading
+  // and the mesh size.
   cout << out.second;
   if (out.second != true)
   {
     warning("This demo is a complex nonlinear problem. Convergence is not guaranteed when modifying some parameters or using PETSC 3.2.");
   }
+
   // Save solution in VTK format
   File file("displacement.pvd");
   file << u;
 
   // Make plot windows interactive
   interactive();
-#endif
+
 #endif
 
  return 0;

demo/undocumented/contact-vi-snes/python/demo_contact-vi-snes.py

-"""This demo program uses of the interface to SNES solver for variational inequalities
- to solve a contact mechanics problems in FEnics.
-The example considers a heavy hyperelastic circle in a box of the same size"""
+"""This demo program uses of the interface to SNES solver for
+ variational inequalities to solve a contact mechanics problems in
+ FEnics.  The example considers a heavy hyperelastic circle in a box
+ of the same size"""
 # Copyright (C) 2012 Corrado Maurini
 #
 # This file is part of DOLFIN.
     print "DOLFIN must be compiled with PETSc to run this demo."
     exit(0)
 
-if not has_cgal():
-    print "DOLFIN must be compiled with CGAL to run this demo."
-    exit(0)
-
 # Create mesh
-mesh = UnitCircleMesh(30)
+mesh = Mesh("../circle_yplane.xml.gz")
 
 V = VectorFunctionSpace(mesh, "Lagrange", 1)
 
 
 # The displacement u must be such that the current configuration x+u
 # remains in the box [xmin,xmax] x [umin,ymax]
-constraint_u = Expression( ("xmax - x[0]","ymax - x[1]"), xmax=1 + DOLFIN_EPS,  ymax=1.0)
-constraint_l = Expression( ("xmin - x[0]","ymin - x[1]"), xmin=-1 - DOLFIN_EPS, ymin=-1.0)
+constraint_u = Expression(("xmax - x[0]","ymax - x[1]"), \
+                          xmax=1.0+DOLFIN_EPS,  ymax=1.0)
+constraint_l = Expression(("xmin - x[0]","ymin - x[1]"), \
+                          xmin=-1.0-DOLFIN_EPS, ymin=-1.0)
 umin = interpolate(constraint_l, V)
 umax = interpolate(constraint_u, V)
 

demo/undocumented/contact-vi-tao/cpp/main.cpp

 // Copyright (C) 2012 Corrado Maurini
-// 
+//
 // This file is part of DOLFIN.
-// 
+//
 // DOLFIN is free software: you can redistribute it and/or modify
 // it under the terms of the GNU Lesser General Public License as published by
 // the Free Software Foundation, either version 3 of the License, or
 // (at your option) any later version.
-// 
+//
 // DOLFIN is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 // GNU Lesser General Public License for more details.
-// 
+//
 // You should have received a copy of the GNU Lesser General Public License
 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
-// 
+//
 // Modified by Corrado Maurini 2013
-// 
+//
 // First added:  2012-09-03
 // Last changed: 2013-04-15
 //
-// This demo program uses of the interface to TAO solver for variational inequalities 
-// to solve a contact mechanics problem in FEnics. 
-// The example considers a heavy elastic circle in a box of the same diameter
+// This demo program uses of the interface to TAO solver for
+// variational inequalities to solve a contact mechanics problem in
+// FEnics.  The example considers a heavy elastic circle in a box of
+// the same diameter
 
 #include <dolfin.h>
 #include "Elasticity.h"
 
 using namespace dolfin;
 
-int main()
+// Lower bound for displacement
+class LowerBound : public Expression
 {
-  #ifdef HAS_TAO
-    
-  // Lower bound for displacement
-  class LowerBound : public Expression
+public:
+  LowerBound() : Expression(2) {}
+  void eval(Array<double>& values, const Array<double>& x) const
   {
-  public:
-    
-    LowerBound() : Expression(2) {}
-    
-    void eval(Array<double>& values, const Array<double>& x) const
-    {
-      double xmin = -1.-DOLFIN_EPS;
-      double ymin = -1.;
-      values[0] = xmin-x[0];
-      values[1] = ymin-x[1];
-    }
-    
-  };
-    
-  // Upper bound for displacement
-  class UpperBound : public Expression
+    const double xmin = -1.0 - DOLFIN_EPS;
+    const double ymin = -1.0;
+    values[0] = xmin - x[0];
+    values[1] = ymin - x[1];
+  }
+};
+
+// Upper bound for displacement
+class UpperBound : public Expression
+{
+public:
+  UpperBound() : Expression(2) {}
+  void eval(Array<double>& values, const Array<double>& x) const
   {
-  public:
-    
-    UpperBound() : Expression(2) {}
-
-    void eval(Array<double>& values, const Array<double>& x) const
-    {
-      double xmax = 1.+DOLFIN_EPS;
-      double ymax = 1.;
-      values[0] = xmax-x[0];
-      values[1] = ymax-x[1];
-    }
-    
-  };
-
-  // Read mesh and create function space
-  UnitCircleMesh mesh(50);
+    const double xmax = 1.0 + DOLFIN_EPS;
+    const double ymax = 1.0;
+    values[0] = xmax-x[0];
+    values[1] = ymax-x[1];
+  }
+};
+
+int main()
+{
+  #ifdef HAS_TAO
+  // Read mesh
+  Mesh mesh("../circle_yplane.xml.gz");
+
+  // Create function space
   Elasticity::FunctionSpace V(mesh);
 
   // Create right-hand side
   Function usol(V);
   Function Ulower(V); // Create a function
   Function Uupper(V); // Create a function
-  
+
   // Assemble the matrix and vector for linear elasticity
   PETScMatrix A;
   PETScVector b;
   assemble(A, a);
   assemble(b, L);
-    
-  // Interpolate expression for Upper bound 
+
+  // Interpolate expression for Upper bound
   UpperBound xu_exp;
   Function xu_f(V);
   xu_f.interpolate(xu_exp);
- 
-  // Interpolate expression for Upper bound 
+
+  // Interpolate expression for Upper bound
   LowerBound xl_exp;
   Function xl_f(V);
-  xl_f.interpolate(xl_exp);  
-    
+  xl_f.interpolate(xl_exp);
+
   // Create the PetscVector associated to the functions
   PETScVector& x  = (*usol.vector()).down_cast<PETScVector>(); // Solution
   PETScVector& xl = (*xl_f.vector()).down_cast<PETScVector>(); // Lower bound
   PETScVector& xu = (*xu_f.vector()).down_cast<PETScVector>(); // Upper bound
 
-  // Solve the problem with the TAO Solver 
+  // Solve the problem with the TAO Solver
   TAOLinearBoundSolver TAOSolver("tao_tron","tfqmr");
-  
+
   // Set some parameters
   TAOSolver.parameters["monitor_convergence"]=true;
   TAOSolver.parameters["report"]=true;
   TAOSolver.parameters("krylov_solver")["monitor_convergence"]=false;
-  
-  // Solve the problem 
+
+  // Solve the problem
   TAOSolver.solve(A, x, b, xl, xu);
-  
-  
+
   // Plot solution
   plot(usol, "Displacement", "displacement");
-  
+
   // Make plot windows interactive
   interactive();
 
   cout << "This demo requires DOLFIN to be configured with TAO" << endl;
 
   #endif
-  
+
  return 0;
 }

demo/undocumented/contact-vi-tao/python/demo_contact-vi-tao.py

     print "DOLFIN must be compiled with TAO to run this demo."
     exit(0)
 
-# Create mesh (use cgal if available)
-mesh = UnitCircleMesh(50)
+# Read mesh
+mesh = Mesh("../circle_yplane.xml.gz")
 
 # Create function space
 V = VectorFunctionSpace(mesh, "Lagrange", 1)

dolfin/nls/PETScSNESSolver.cpp

       }
       else
       {
-        #if PETSC_HAVE_MUMPS
-        lu_method = "mumps";
+        #if PETSC_HAVE_SUPERLU_DIST
+        lu_method = "superlu_dist";
         #elif PETSC_HAVE_PASTIX
         lu_method = "pastix";
-        #elif PETSC_HAVE_SUPERLU_DIST
-        lu_method = "superlu_dist";
+        #elif PETSC_HAVE_MUMPS
+        lu_method = "mumps";
         #else
         dolfin_error("PETScSNESSolver.cpp",
                      "solve linear system using PETSc LU solver",
     {
       // Here, x is the model vector from which we make our Vecs that
       // tell PETSc the bounds.
-      Vec ub;
-      Vec lb;
+      Vec ub, lb;
 
       PETScVector dx = x.down_cast<PETScVector>();
       VecDuplicate(*dx.vec(), &ub);