Mixed Dirichlet and Neumann boundary conditions

[bjpal...@gmail.com]

Hi,

I'm an FE newbie and I'm trying to solve an electrostatic problem with a spatially varying dielectric function and no charges (div epsilon grad phi = 0). The field is determined by specifying a combination of Dirichlet and Neumann conditions on the boundary. The boundary is divided into disjoint segments with Neumann conditions (vanishing normal gradient) on some parts of the boundary and constant values (Dirichlet conditions) on the other disjoint segments of the boundary. I'm trying to set the Dirichlet conditions based on boundary IDs and then pass these values into the matrix equations using the ConstraintMatrix object. The boundary IDs run from 0-3 with 3 representing the Neumann condition and 0-2 representing different Dirichlet values. At the moment I'm just using a regular grid so there is no issue with hanging nodes or anything like that. The code I've been using looks like

      p_constraints.clear();
      typename FunctionMap<3>::type dirichlet_boundary_functions;
      ConstantFunction<3> gate(data.gatePotential(),1);
      dirichlet_boundary_functions[0] = &gate;
      ConstantFunction<3> source(data.sourcePotential(),1);
      dirichlet_boundary_functions[1] = &source;
      ConstantFunction<3> drain(data.drainPotential(),1);
      dirichlet_boundary_functions[2] = &drain;

      VectorTools::interpolate_boundary_values(p_dof,
          0, gate, p_constraints);
      VectorTools::interpolate_boundary_values(p_dof,
          1, source, p_constraints);
      VectorTools::interpolate_boundary_values(p_dof,
          2, drain, p_constraints);

      p_constraints.close();

where p_constraints is an object of type ConstraintMatrix. I've used p_constraints to create a sparsity pattern and then used that to create the system matrix. I've also used p_constraints to modify the solution vector using

      p_constraints.distribute(p_solution);

and tried to add the constraints to the system matrix and right hand side using the call

        p_constraints.distribute_local_to_global(cell_matrix, cell_rhs,
            local_dof_indices,
            p_system, p_rhs);

when looping over FE cells. However, I don't think anything is happening. As far as I can tell, the surface terms in the weak formulation just drop out (I'm not using a coupled potential-displacement field formulation) so if the constraints don't show up in the equations, you just end up with the equation Ax = 0, which isn't very interesting. What do I need to do to get the constraints into my equations?

Bruce Palmer