assemble matrix after mesh refinement

[Mauro Murer]

Hello everyone!

I'm trying to rewrite Step-35 (Navier-Stokes tutorial) without splitting the implementation for component (In Step 35 one thread resolve the x component and the other one the y component).

Proceeding in this way I can use "AffineConstraints" class to manipulate constraints in the system matrices.

All went well until I tried to implement a mesh refinement.

In particular, when I assemble the gradient operator after the refinement of the mesh (increasing the dof of the problem), this error occurs when "distribute_local_to_global" function is called:

The violated condition was:
    (this_cols[counter] == col_indices[i]) || (values[i] == number2())
Additional information:
    You are trying to access the matrix entry with index <35518,4519>, but this entry does not exist in the sparsity pattern of this matrix.

The most common cause for this problem is that you used a method to build the sparsity pattern that did not (completely) take into account all of the entries you will later try to write into. An example would be building a sparsity pattern that does not include the entries you will write into due to constraints on degrees of freedom such as hanging nodes or periodic boundary conditions. In such cases, building the sparsity pattern will succeed, but you will get errors such as the current one at one point or other when trying to write into the entries of the matrix.

I found in the documentation the "add_entries" command, in such a way to add the matrix entry in the sparsity pattern, but I don't know if it is the right way to proceed. Any suggestions???

Thanks a lot

The piece of the code, where the gradient operator is assembled, is reported below:

    /* Gradient operator*/

    DynamicSparsityPattern dsp3(dof_handler_velocity.n_dofs(),
                                dof_handler_pressure.n_dofs());
    DoFTools::make_sparsity_pattern(dof_handler_velocity,
                                    dof_handler_pressure,
                                    dsp3);
    sparsity_pattern_pres_vel.copy_from(dsp3);
    pres_Diff.reinit(sparsity_pattern_pres_vel);

    unsigned int vel_dpc(fe_velocity.dofs_per_cell);
    unsigned int pres_dpc(fe_pressure.dofs_per_cell);

    FullMatrix<double> local_grad(fe_velocity.dofs_per_cell, fe_pressure.dofs_per_cell);

    std::vector<types::global_dof_index> vel_local_dof_indices(fe_velocity.dofs_per_cell);
    std::vector<types::global_dof_index> pres_local_dof_indices(fe_pressure.dofs_per_cell);

    unsigned int nqp(quadrature_velocity.size());

    FEValues<dim> fe_val_vel(fe_velocity, quadrature_velocity,
                             update_values | update_gradients | update_JxW_values);

    FEValues<dim> fe_val_pres(fe_pressure, quadrature_velocity,
                              update_values);

    auto vel_cell = dof_handler_velocity.begin_active();
    const auto endc = dof_handler_velocity.end();
    auto pres_cell = dof_handler_pressure.begin_active();

    for (; vel_cell != endc; ++vel_cell, ++pres_cell) {

        fe_val_vel.reinit(vel_cell);
        fe_val_pres.reinit(pres_cell);

        local_grad = 0;

        for (unsigned int q_point = 0; q_point < nqp; ++q_point) {
            for (unsigned int i = 0; i < vel_dpc; ++i) {
                const unsigned int component_i = fe_velocity.system_to_component_index(i).first;
                for (unsigned int j = 0; j < pres_dpc; ++j) {

                    local_grad(i, j) +=
                            fe_val_pres.shape_value(j, q_point)*
                            fe_val_vel.shape_grad(i, q_point)[component_i] *
                            (-fe_val_vel.JxW(q_point));
                }
            }
        }

        vel_cell->get_dof_indices(vel_local_dof_indices);
        pres_cell->get_dof_indices(pres_local_dof_indices);

        velocity_hanging_node_constraints.distribute_local_to_global(local_grad, vel_local_dof_indices,
                                                                     pressure_hanging_node_constraints,
                                                                     pres_local_dof_indices, pres_Diff);


    }

[Wolfgang Bangerth]

On 10/15/19 1:55 AM, Mauro Murer wrote:
> s in the system matrices.
> All went well until I tried to implement a mesh refinement.
> In particular, when I assemble the gradient operator after the refinement of
> the mesh (increasing the dof of the problem), this error occurs when
> "distribute_local_to_global" function is called:
>

> *The violated condition was:

>     (this_cols[counter] == col_indices[i]) || (values[i] == number2())
> Additional information:
> You are trying to access the matrix entry with index <35518,4519>, but this
> entry does not exist in the sparsity pattern of this matrix.
>
> The most common cause for this problem is that you used a method to build the
> sparsity pattern that did not (completely) take into account all of the
> entries you will later try to write into. An example would be building a
> sparsity pattern that does not include the entries you will write into due to
> constraints on degrees of freedom such as hanging nodes or periodic boundary
> conditions. In such cases, building the sparsity pattern will succeed, but you
> will get errors such as the current one at one point or other when trying to

> write into the entries of the matrix.*

>
>
> I found in the documentation the "add_entries" command, in such a way to add
> the matrix entry in the sparsity pattern, but I don't know if it is the right
> way to proceed. Any suggestions???

You need to take constraints into account when you build the sparsity pattern
but I believe from your screen shot that you don't. Take a look at step-6, for
example:
https://github.com/dealii/dealii/blob/master/examples/step-6/step-6.cc#L225

Best
W.

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

[Mauro Murer]

Thanks a lot for your reply.

I tried to take constraints into account in the sparsity pattern, but I don't know how to do it for a rectangular matrix such in my case.

I have rows equal to the number of velocity dofs and columns equal to pressure dofs.

I didn't found any function for sparsity pattern in the documentation that take in account different constraints (in my case velocity and pressure constraints).

Maybe I didn't looking for enough, but I actually stuck in this problem.

Thanks

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3efefd13-6554-9fb3-fd7e-b6bf230d4b00%40colostate.edu.

[Wolfgang Bangerth]

On 10/21/19 2:55 AM, Mauro Murer wrote:
>
> I tried to take constraints into account in the sparsity pattern, but I don't
> know how to do it for a rectangular matrix such in my case.
> I have rows equal to the number of velocity dofs and columns equal to pressure
> dofs.
> I didn't found any function for sparsity pattern in the documentation that
> take in account different constraints (in my case velocity and pressure
> constraints).

Then just build the sparsity pattern as you do and use
AffineConstraints::compress() after the fact -- I believe that there is such a
function that can work on rectangular matrices.

In general, you are discussing why it is simpler to just keep all components
of the solution in the same DoFHandler object. In that case, the rectangular
matrix is just one block in a square matrix, and all of the constraint
handling things work as expected.

[Mauro Murer]

Actually I can't find the function AffineConstraints::compress().

I tried to use SparseMatrix::compress(::VectorOperation::values), on the other hand, but the error still remain.

Maybe I don't use the right function??

Thanks

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/146c2446-0d06-1cec-a908-0b474e4747dd%40colostate.edu.

[Wolfgang Bangerth]

On 10/21/19 8:32 AM, Mauro Murer wrote:
> Actually I can't find the function AffineConstraints::compress().

Sorry, I meant ::condense().

[Mauro Murer]

It seems it doesn't work for rectangular objects!!!

An error occurred in line <1191> of file </home/murer/Desktop/dealii/dealii-9.1.1/include/deal.II/lac/affine_constraints.templates.h> in function

    void dealii::AffineConstraints<number>::condense(dealii::DynamicSparsityPattern&) const [with number = double]
The violated condition was:
    sparsity.n_rows() == sparsity.n_cols()
Additional information:
    This function only works for quadratic objects!

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4d17534b-c63e-5d49-d86e-0151b3ee0dfe%40colostate.edu.

[Wolfgang Bangerth]

Mauro,
apologies for not getting back to this earlier. Have you found a solution?

I think what you find is something that we had to learn over time and have
learned since step-35 was written: That it's typically easiest to just put all
solution variables into one large finite element system, and then build a
square matrix that, if decomposed, may have rectangular matrices as building
blocks. If we were to re-write step-35 today, that's the design we'd choose --
in fact, if anyone wanted to re-write the program in this style, I think that
would make for an excellent project we would very gladly take!

So that would be what I suggest you consider: Rewriting things in such a way
that the various matrices become blocks in a large square matrix. If you want
to go the cheap way, you can always just build a large square matrix for which
you really only care about the one block in question and the rest is only
there to make sure that you have a square matrix, rather than containing
anything useful.

Best
W.

On 10/21/19 9:22 AM, Mauro Murer wrote:
> It seems it doesn't work for rectangular objects!!!
>
> An error occurred in line <1191> of file
> </home/murer/Desktop/dealii/dealii-9.1.1/include/deal.II/lac/affine_constraints.templates.h>
> in function
>     void
> dealii::AffineConstraints<number>::condense(dealii::DynamicSparsityPattern&)
> const [with number = double]
> The violated condition was:
>     sparsity.n_rows() == sparsity.n_cols()
> Additional information:
>     This function only works for quadratic objects!
>
> Il giorno lun 21 ott 2019 alle ore 16:35 Wolfgang Bangerth

> <bang...@colostate.edu <mailto:bang...@colostate.edu>> ha scritto:

>
> On 10/21/19 8:32 AM, Mauro Murer wrote:
> > Actually I can't find the function AffineConstraints::compress().
>
> Sorry, I meant ::condense().
>
> Best
>   W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth          email: bang...@colostate.edu

> <mailto:bang...@colostate.edu>

>                             www: http://www.math.colostate.edu/~bangerth/
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com

> <mailto:dealii%2Bunsu...@googlegroups.com>.

> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/4d17534b-c63e-5d49-d86e-0151b3ee0dfe%40colostate.edu.
>

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com

> <mailto:dealii+un...@googlegroups.com>.

> To view this discussion on the web visit

> https://groups.google.com/d/msgid/dealii/CAK3LiWecHAVGMq%2Bj3GXwwJ%3DJqat-1OY4PTFWZgH_X0FFA5EbbA%40mail.gmail.com
> <https://groups.google.com/d/msgid/dealii/CAK3LiWecHAVGMq%2Bj3GXwwJ%3DJqat-1OY4PTFWZgH_X0FFA5EbbA%40mail.gmail.com?utm_medium=email&utm_source=footer>.