Which b.c.s to reproduce 1-D solution with 2-D code

There is a subtly here related to the Neumann boundary condition. DOLFIN interpolates the BCs in finite element spaces, which can lead to some 'creep' around corners if the domain is not properly specified. What you need to do is mark the boundary on which you're applying the force bc, e.g

In the UFL file, do something like:

  L = f*vdx + h*v*ds(3)

This means the h*v should only be integrated over the boundary subdomain '3'.

In the C++ code:

  // Right boundary
  class RightBoundary : public SubDomain
  {
    bool inside(const Array<double>& x, bool on_boundary) const
    {
      if (1.0 - x[0] < DOLFIN_EPS && on_boundary)
        return true;
      else
        return false;
    }
  };

  // Create subdomain
  RightBoundary right_boundary;

  // Create mesh function over facets (FacetFunction is new in 0.9.9, otherwise use
  // MeshFunction<unsigned int> right_boundary_function(mesh, 1) for a 2D problem);
  FacetFunction<unsigned int> right_boundary_function(mesh);

  // Mark the boundary with '3'
  right_boundary.mark(right_boundary_function, 3);

If you're using the VariationalProblem class, then do

  // Create variational problem
  VariationalProblem pde(a_form, L, bc, 0, &right_boundary_function, 0);

This will pass the subdomains to the assembler. If you're calling the assembler directly, take a look in

  dolfin/fem/Assembler.h

for how to pass subdomains.

The elastodynamics demo also shows how to use boundary subdomains when applying force bcs (both in C++ and Python).