non-homogeneous dirichlet boundary with constraint matrix. Get zero solution on the boundary.

[Anders Ström]

Hello,

I'm trying to solve a 2x2 block system with non-homogenous Dirichlet boundary condition for the first component of the solution. The solution I get is zero on the boundary though. The b.c. is set with a constraint matrix and a component mask like this, in make_grid_and_dof,

    constraints.clear ();
    FEValuesExtractors::Scalar ys(0);
    ComponentMask y_mask = fe.component_mask(ys);
    DoFTools::make_hanging_node_constraints (dof_handler,
                                             constraints);
    VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
                                              BoundaryValues<dim>() ,//ZeroFunction<dim>(2)
                                              constraints,
                                              y_mask);
    constraints.close ();

and the boundary function class looks like this,

class BoundaryValues : public Function<dim>
{
public:
    BoundaryValues() : Function<dim>(2) {}
    virtual double value (const Point<dim>  &p,
                          const unsigned int component = 0) const;
    virtual void vector_value (const Point<dim> &p,
                                       Vector<double>   &values) const;
};

template<int dim>
void
BoundaryValues<dim>::vector_value (const Point<dim> &p,
                                   Vector<double>   &values) const
{
  for (unsigned int c=0; c<this->n_components; ++c)
    values(c) = BoundaryValues<dim>::value (p, c);
}

template<int dim>
double BoundaryValues<dim>::value(const Point<dim>  &p,
                                  const unsigned int component) const
{
    double x=p[0], y=p[1];
    double return_value;
    return_value=0;
    if(component == 0)
        if(x<0.5 && y<0.5){
            return_value = (2*x-1)*(2*x-1)*(2*y-1)*(2*y-1);
            std::cout<< component << ", " << x << ", " << y << "," << return_value<<std::endl;
        }
    return return_value;
}

In assemble_system the I then use,

cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
                                        local_dof_indices,
                                        system_matrix,
                                        system_rhs);

The system is a PDE-constrained optimization problem with Poisson equation. block(0) of the solution corresponds to the state for which the boundary condition should be set. block(1) is the control and should not have any boundary conditions.

I have attached the code if there is no apparent issue in the above. Thankful to any answer.

Regards,
Anders  








[Anders Ström]

You should call,

constraints.distribute(solution);

after solving the system.

[Jean-Paul Pelteret]

Just to confirm, is this the solution to your problem? It seems likely that this is what you need to do.

[Anders Ström]

yes, this solved it. I was going to marked as solved or close it but I couldn't figure out how.

[Jean-Paul Pelteret]

No problem! I'm glad that you managed to figure it out :-)