Imposing homogeneous Dirichlet boundary conditions with apply_boundary_values

[Ali Seddiq]

Dear all,

I have a problem with imposing homogeneous Dirichlet boundary conditions for pressure in a vector-valued problem.

I have defined interpolate_ and apply_boundary_values as following  in the assembly routine:

std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,

                                                                     1,

                                                                     BoundaryValues<dim>(), //copied from step-22 with slight modifications
                                                                     boundary_values,
                                                                    fe.component_mask(pressure));

  MatrixTools::apply_boundary_values(boundary_values, system_matrix, solution, system_rhs);

But it throws an exception (in line 261, matrix_tools.cc) stating:

 matrix.get_sparsity_pattern().get_column_indices() == solution.get_block_indices()
Additional information:
    You are providing a matrix whose subdivision into blocks in either row or column direction does not use the same block sizes as the solution vector or right hand side vectors, respectively.

(Stacktrace:...
 void dealii::MatrixTools::apply_boundary_values<double>(std::map<unsigned int, double, std::less<unsigned int>, std::allocator<std::pair<unsigned int const, double> > > const&, dealii::BlockSparseMatrix<double>&, dealii::BlockVector<double>&, dealii::BlockVector<double>&, bool)

)

The code (and setup_system in particular) is basically based on step-20 (i.e. I used the component_wise renumbering approach, BlockSparseMatrix and BlockVectors, but no Schur complement, no LinearOperator, and a matrix much simpler). But I still do not understand the problem with my definitions. Could anyone have some hints please?

Thank you very much,

Ali

[Wolfgang Bangerth]

Ali,

> I have a problem with imposing homogeneous Dirichlet boundary conditions for
> pressure in a vector-valued problem.
>
> I have defined interpolate_ and apply_boundary_values as following  in the
> assembly routine:
>
> std::map<types::global_dof_index, double> boundary_values;
> VectorTools::interpolate_boundary_values(dof_handler,
>                                                                      1,
>
> BoundaryValues<dim>(), //copied from step-22 with slight modifications
>
> boundary_values,
>
> fe.component_mask(pressure));
>
>   MatrixTools::apply_boundary_values(boundary_values, system_matrix,
> solution, system_rhs);
>
> But it throws an exception (in line 261, matrix_tools.cc) stating:
>
>  matrix.get_sparsity_pattern().get_column_indices() ==
> solution.get_block_indices()
> Additional information:
>     You are providing a matrix whose subdivision into blocks in either row or
> column direction does not use the same block sizes as the solution vector or
> right hand side vectors, respectively.

Well, what is the subdivision of matrix and vector, and why are they
different? How did you subdivide them?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Ali Seddiq]

Dear Wolfgang,

Thank you for your reply. In order to reconstruct the following final discretized form as

(M+A)u + A'p=F+G,

 a 2*2 matrix is subdivided into blocks in a way that (block)components (0,0) and (1,1) are M+A and A' respectively and components  (0,1) and (1,0) set to zero.

Solution vector is declared as a BlockVector but contains/subdivided into two scalar components u and p.  Right hand side is also a BlockVector containing F_u and G_p.

I hope I could have addressed your point here. Would you please advise me if there are any discrepancies in my structuring?

Thank you very much,

Ali

[Wolfgang Bangerth]

On 2/25/22 19:44, Ali Seddiq wrote:
>
> Thank you for your reply. In order to reconstruct the following final
> discretized form as
> (M+A)u + A'p=F+G,
>  a 2*2 matrix is subdivided into blocks in a way that (block)components (0,0)
> and (1,1) are M+A and A' respectively and components  (0,1) and (1,0) set to zero.
> Solution vector is declared as a BlockVector but contains/subdivided into two
> scalar components u and p.  Right hand side is also a BlockVector containing
> F_u and G_p.
> I hope I could have addressed your point here. Would you please advise me if
> there are any discrepancies in my structuring?

I have no idea. The error message says

You are providing a matrix whose subdivision into blocks in either row or
column direction does not use the same block sizes as the solution vector or
right hand side vectors, respectively.

You need to find out why that is so. Some debug output right before you call
MatrixTools::apply_boundary_values() will be able to tell you whether the
error message is correct or not. If the error message is correct, you ought to
find out why the subdivision is how it is.