Mixed codimensionality

[Alex Quinlan]

Dear all,

I've been working with deal.ii in the solid mechanics field for both 3D and 2D cases.

I am now being asked by my colleagues to develop a model that can combine elements of different dimensionality.  This is approach is already implemented in my colleagues' abaqus models, and so we're looking to replicate the functionality in deal.ii.

The elements to be combined would be FE_System elements of type:

- 3D volumetric element in 3D space <dim=3, spacedim=3>

- 2D planar element in 3D space <dim=2, spacedim=3>

- 1D linear element in 3D space <dim=1, spacedim=3>

I've attached a sketch of an application that might use this. 

From what I've seen, deal.ii does not like to mix <dim>'s.  Does anyone have thoughts on how I might go about implementing a model with these elements?  If I need to modify or implement new functionalities, could you estimate the magnitude of time required to make changes? (i.e. days, weeks, months)

I realize that there are alternative methods to represent the physical system that could avoid this problem.  However, I am essentially being asked to automate the conversion from existing abaqus input files into deal.ii.  I appreciate any ideas you all might have.

Thanks,

Alex

[Daniel Arndt]

Alex,

FiniteElements in deal.II are used with a Triangulation and DoFHandler of the same dimension and space dimension. 

From your sketch, it seems like the Triangulation for your 3 types of finite elements is (expectedly) different. Thus, it doesn't make much sense to combine those in a FESystem.

It's probably better to define three separate Triangulations and use the respective finite element classes for them. The tricky would likely be getting the constraints between them right.

What does your weak form look like anyway?

Best,

Daniel

--
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/a7ff55e9-b406-403b-b176-75690c2d121en%40googlegroups.com.

[Alex Quinlan]

Hi Daniel,

Thanks for the response.  I am essentially using the approach in Step-8, where I am seeking a vector-valued displacement solution using FE_System.  So, I've been using:

dealii::FESystem<3>(FE_Q<3>(1), 3 )

and

dealii::FESystem<2,3>(FE_Q<2,3>(1), 3 )

Based on your suggestion, I am picturing doing something like this:

// Shared Objects

AffineConstraints<double> constraints;
SparsityPattern    sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;    // This is the solution vector
Vector<double> system_rhs;  // This is the load / right-hand-side vector

// Dimension dependent objects

// 3D

Triangulation<3,3>   triangulation_3D;
DoFHandler<3,3>      dof_handler_3D;

FESystem<3,3>        fe_3D;

// 2D

Triangulation<2,3>   triangulation_2D;
DoFHandler<2,3>      dof_handler_2D;

FESystem<2,3>        fe_2D;

// 1D

Triangulation<1,3>   triangulation_1D;
DoFHandler<1,3>      dof_handler_1D;

FESystem<1,3>        fe_1D;

Then I need to proceed with each of these triangulations, and somehow link them together within the system_matrix and system_rhs.

I have been building the system matrix using this line for each cell in the dof_handler.active_cell_iterators:

constraints.distribute_local_to_global(      cell_matrix, cell_rhs, local_dof_indices,      
                                             system_matrix, system_rhs);

I see two challenges using this:

1)  keeping track of the dof_indices across 3 dof_handlers so they can be added to the correct place in the system_matrix

2) Connecting the triangulations within the system matrix.

Have I understood your suggestion correctly?  Do you forsee any other challenges?

Thanks very much,

Alex

[Daniel Arndt]

Alex,

If you can formulate your problem/approach in a way (iteratively?) that you could solve separately on the three triangulations (using input from the solution on the other triangulations), that might be easier.

On the other hand, if you only want one system matrix you would identify matching dofs manually and merge the system matrices manually.

Maybe, someone else has better ideas and more experience than me. I'm still curious how your full weak form looks like.

Best,

Daniel

--
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/7b15f86c-0069-42c9-b7ed-4bde752ee07bn%40googlegroups.com.

[Wolfgang Bangerth]

On 10/27/23 12:21, Daniel Arndt wrote:
>
> If you can formulate your problem/approach in a way (iteratively?) that you
> could solve separately on the three triangulations (using input from the
> solution on the other triangulations), that might be easier.

That corresponds to an operator splitting scheme or, in the language of linear
solvers, an Uzawa iteration.

The alternative is to build a block system in which each block corresponds to
one triangulation (2 or 3 dimensional). You'd build each block in the same way
as you would with a single triangulation, including using the same indexing of
rows/columns that corresponds to the DoFHandlers on these triangulations. The
tricky part is assembling the off-diagonal blocks where you'd need to work a
bit harder. If you can assemble such a linear system, then the solution
becomes both relatively easy and efficient because you can use a coupled
solver scheme.

You might also want to look at some of the "non-matching" tutorial programs to
see inspirations for problems that use different triangulations at once.

Best
W.

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

[Alex Quinlan]

Dear Wolfgang and Daniel,

Thanks for your directions.  I've reviewed Step-20, Step-60, and Step-70 and started working on my own example.
In this case, I'm trying to combine two types of dimensionality: the <3,3> Solid Elements and the <2,3> "Shell" Elements. I've attached a picture of the meshes.  I am trying to do a block system.  For the moment, I am not trying to link the two meshes (off-diagonal blocks), only get them into the same system matrix. Each is independently constrained with appropriate boundary conditions.  I've run the meshes separately and they work fine.

My understanding is that I cannot combine several objects because of the difference in <dim,spacedim>.  So, I've created duplicates of several objects, using "_sld" to note for the Solid and "_sh" for the shell.

Triangulation<3> triangulation_sld;  
Triangulation<2,3> triangulation_sh;
DoFHandler<3>    dof_handler_sld;
DoFHandler<2,3>    dof_handler_sh;    
FESystem<3>      fe_sld;
FESystem<2,3>      fe_sh;
const MappingFE<3>    mapping_sld;
const MappingFE<2,3>  mapping_sh;
const QGauss<3> quadrature_sld;
const QGauss<2> quadrature_sh;

I then create the following objects that will be shared by both the Solid and Shell elements.

AffineConstraints<double> constraints;    
BlockSparsityPattern      sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> preconditioner_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
BlockVector<double> locally_relevant_solution;

I proceed by reading two separate external meshes and attaching them to the appropriate triangulations.
Then I begin the Setup_System().  I start by getting the number of DOFs from each of the triangulations and adding them
in a vector of the DOFs for each block.  I am saying that the Solid is block 0 and the Shell is block 1.

const unsigned int n_dof_sld = dof_handler_sld.n_dofs();  // = 24
const unsigned int n_dof_sh = dof_handler_sh.n_dofs();    // = 48
const std::vector<unsigned int> dofs_per_block = {n_dof_sld, n_dof_sh};

I then begin work on the sparsity pattern:

BlockDynamicSparsityPattern dsp(dofs_per_block , dofs_per_block);
for (unsigned int i=0; i<dofs_per_block.size(); ++i)
    for (unsigned int j=0; j<dofs_per_block.size(); ++j)
    dsp.block(i,j).reinit(dofs_per_block[i], dofs_per_block[j]);

dsp.collect_sizes();

sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
system_rhs.reinit(dofs_per_block);

From what I can see, the system_matrix looks right with n_blocks = 2 and  start_indices = {0, 24, 72}, as expected.

Next, I run through the Assemble_System() to create a local cell matrix. I start with the Solid elements and cell_matrix_sld. I'm skipping the calculation of the cell matrix for brevity's sake.

for (const auto &cell : dof_handler_sld.active_cell_iterators())
    {
    // *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix ***

    cell->get_dof_indices(local_dof_indices_sld);
    // Add the cell matrix to the system matrix
    for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
          for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)
            system_matrix.add(local_dof_indices_sld[i],
                              local_dof_indices_sld[j],
                              cell_matrix_sld(i, j));
    }  // end of cell loop

The bit of code 'system_matrix.add(...)' works for when i=0 and j=0, but fails on the next loop when j=1.  I get this error:


An error occurred in line <1931> of file <.../deal.II/lac/sparse_matrix.h> in function
    void dealii::SparseMatrix<Number>::add(dealii::SparseMatrix<Number>::size_type, dealii::SparseMatrix<Number>::size_type, number) [with number = double; dealii::SparseMatrix<Number>::size_type = unsigned int]
The violated condition was:
    (index != SparsityPattern::invalid_entry) || (value == number())
Additional information:
    You are trying to access the matrix entry with index <0,1>, 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.

So, it appears that I'm not setting something up correctly with the system_matrix and sparsity pattern.  Step-20 uses the same system_matrix.add(...) approach, and does not have any additional preparation of the system_matrix that I can see.

One thing I skipped is this line from Step 20:

DoFTools::make_sparsity_pattern(dof_handler, dsp);

I'm unable to use this line because I have the two separate DOF handlers.  Do you have any idea of what I can do here?

Best regards,

Alex

[Wolfgang Bangerth]

Alex:

> [...]

I think everything up to this point looks reasonable.

> From what I can see, the system_matrix looks right with n_blocks = 2 and
>  start_indices = {0, 24, 72}, as expected.
> Next, I run through the Assemble_System() to create a local cell matrix. I
> start with the Solid elements and cell_matrix_sld. I'm skipping the
> calculation of the cell matrix for brevity's sake.
>
> for (const auto &cell : dof_handler_sld.active_cell_iterators())
>     {
>     // *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix ***
>     cell->get_dof_indices(local_dof_indices_sld);
>     // Add the cell matrix to the system matrix
>     for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
>           for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)
>             system_matrix.add(local_dof_indices_sld[i],
>                               local_dof_indices_sld[j],
>                               cell_matrix_sld(i, j));
>     }  // end of cell loop
>
>
> The bit of code 'system_matrix.add(...)' works for when i=0 and j=0, but fails
> on the next loop when j=1.  I get this error:

The issue here is that you're indexing into the *global* matrix with the
indices of only one of the DoFHandler objects. That can't go well. It *should*
work if you do the indexing only in the matrix block that corresponds to
dof_handler_sld:

system_matrix.block(block_index_sld, block_index_sld)

.add(local_dof_indices_sld[i],
local_dof_indices_sld[j],
cell_matrix_sld(i, j));

where you define block_index_sld somewhere to either be zero or one.

In your case, this won't work yet because you hadn't built the sparsity
pattern yet. You can do that the same way:

DoFTools::make_sparsity_pattern(dof_handler_sld,
block_dsp.block(block_index_sld, block_index_sld));

[Luca Heltai]

Dear Alex,

if you send me your github user, I’ll give you access to a code that does exactly what you are trying to do, so that you can get some inspiration on how to assemble the coupling matrices.

There is an open PR (the review has been stalling for some time on that) https://github.com/dealii/dealii/pull/15773 that tackle exactly this use case.

Let me know if you need some assistance.

L.

> --
> 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/9a890d55-91f6-ec16-98be-acdb730919d3%40colostate.edu.

[Alex Quinlan]

Dear Luca and Wolfgang,

Thanks for your pointers; I will work to implement them.

My github account is: https://github.com/AlexQuinlan

Thanks,

Alex

[Alex Quinlan]

Alright, I've had some decent success, though I've hit another snag in the solver section.  Your suggestions have opened the doors and actually allowed me to shift gears to my preferred approach.

I've changed the end of Setup_System() to be: 

    dsp.collect_sizes();
    DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0));
    DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1));

    sparsity_pattern.copy_from(dsp);
    system_matrix.reinit(sparsity_pattern);
    system_rhs.reinit(dofs_per_block);
    solution.reinit(dofs_per_block);
    std::cout << " " << std::endl;

Doing so allowed me to successfully run your suggested code:

        for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
          for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)

            system_matrix.block(0,0).add(local_dof_indices_sld[i],

                              local_dof_indices_sld[j],
                              cell_matrix_sld(i, j));

This success prompted me to try a different approach to applying the cell stiffness matrix to the global stiffness matrix.  I have an object of AffineConstraints, so I wanted to use the AffineConstraints::distribute_local_to_global().

For the block system, I prepared the constraints like this:   (Note: nodal_bcs() is a routine that I've written to apply constraints that I've read in from a .csv file)

   

    // Set up Solid Constraints

    constraints_sld.clear();
    nodal_bcs<3,3>(solidbcfile,  dof_handler_sld   ,constraints_sld);
    constraints_sld.close();

    // Set-up Shell Constraints
    constraints_sh.clear();
    nodal_bcs<2,3>(shellbcfile,  dof_handler_sh   ,constraints_sh);
    constraints_sh.close();

    // Combine Solid and Shell constraints

    constraints.copy_from(constraints_sld);
    constraints_sh.shift(n_dof_sld);
    constraints.merge(constraints_sh );

This allowed me to add this line to the Cell loop in the Assemble_System():

constraints.distribute_local_to_global(
                  cell_matrix_sld, cell_rhs_sld, local_dof_indices_sld, system_matrix.block(0,0), system_rhs.block(0));

All of this now seems to work and the matrices seem to be correct.

-----

Now, I'm having trouble with the solver. I think the cause is the empty off-diagonal matrix  I am using the approach outlined in Step-20:

       const auto &M = system_matrix.block(0, 0);
    const auto &B = system_matrix.block(0, 1);
    const auto &F = system_rhs.block(0);
    const auto &G = system_rhs.block(1);
    auto &U = solution.block(0);
    auto &P = solution.block(1);
    const auto op_M = linear_operator(M);
    const auto op_B = linear_operator(B);

    ReductionControl         reduction_control_M(2000, 1.0e-18, 1.0e-10);
    SolverCG<Vector<double>> solver_M(reduction_control_M);
    PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
    preconditioner_M.initialize(M);

    const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
    const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
    const auto op_aS = transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;

    IterationNumberControl   iteration_number_control_aS(30, 1.e-18);
    SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
    const auto preconditioner_S = inverse_operator(op_aS, solver_aS, PreconditionIdentity());
    const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
    SolverControl            solver_control_S(2000, 1.e-12);
    SolverCG<Vector<double>> solver_S(solver_control_S);
    const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);

    P = op_S_inv * schur_rhs;
    std::cout << solver_control_S.last_step()
              << " CG Schur complement iterations to obtain convergence."
              << std::endl;
    U = op_M_inv * (F - op_B * P);

At the red line above, I am given this error:

An error occurred in line <637> of file </home/.../lac/solver_cg.h> in function
    void dealii::internal::SolverCG::IterationWorker<VectorType, MatrixType, PreconditionerType, <template-parameter-1-4> >::do_iteration(unsigned int) [with VectorType = dealii::Vector<double>; MatrixType = dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload>; PreconditionerType = dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload>; <template-parameter-1-4> = int]
The violated condition was:
    std::abs(p_dot_A_dot_p) != 0.
Additional information:
    A piece of code is attempting a division by zero. This is likely going
    to lead to results that make no sense.

I'm guessing this is related to the inverse operators?  So, I'm curious if there's some way to use this approach for problems that have a Schur complements that can be either zero or non-zero.

Thanks,

Alex

[Alex Quinlan]

Sorry, please disregard my last question. I've reviewed the Step-20 solver that I tried to use and I see that it is a different set-up to my problem.  I'm now working to apply the Schur complement in the correct manner.

[Alex Quinlan]

Dear all,

 

I've come up against another issue in my mixed-dimensional mesh.  I am now trying to implement MPI, which requires partitioning the triangulation. 

In my case, I have two triangulations: one for the solid <3,3> elements and one for the shell <2,3> elements.  These two triangulations are couple at multiple points.  I combine them into a blockmatrix with the form:

[ A_solid   ,    B_coupling;

  C_coupling,    D_shell      ]

I may also need to expand this in the future to a matrix with 3x3 or nxn blocks.

Do you have any thoughts on how I can partition this problem?  The suite of GridTools::partition_triangulation seem like they will not work for my two-triangulation approach.

Best regards,

Alex

[Wolfgang Bangerth]

On 1/16/24 07:03, Alex Quinlan wrote:
>
> I've come up against another issue in my mixed-dimensional mesh.  I am now
> trying to implement MPI, which requires partitioning the triangulation.
>
> In my case, I have two triangulations: one for the solid <3,3> elements and
> one for the shell <2,3> elements.  These two triangulations are couple at
> multiple points.  I combine them into a blockmatrix with the form:
>
> [ A_solid   ,    B_coupling;
>   C_coupling,    D_shell      ]
>
> I may also need to expand this in the future to a matrix with 3x3 or nxn blocks.
>
> Do you have any thoughts on how I can partition this problem?  The suite of
> GridTools::partition_triangulation seem like they will not work for my
> two-triangulation approach.

Can you elaborate what the problem you are facing is? At least for the
diagonal blocks of the matrix, you would simply use a matrix class such as
PETScWrappers::MPI::SparseMatrix that knows how to partition itself onto
different processes. Similarly, for your vectors, you'd use the appropriate
block vector class in which each vector block is partitioned. You'd have to
initialize these individually, but I suspect you're already doing that anyway.
In each case, blocking (into the 2x2 structure you show above) and
partitioning (among MPI processes) are orthogonal issues.

The biggest issue I see is how you can build the off-diagonal blocks above
because to do so, you need to know information about both meshes on the same
process, and generally the partition of one triangulation owned by a specific
process does not geometrically overlap with the partition of the other
triangulation owned by the same process.

[Alex Quinlan]

Hi Wolfgang,

Thanks for your feedback.  I have this two-triangulation problem working in serial, and I have other jobs that work with MPI.  When I put them together, I end up with something like the code below.  It seems like there will be some issues, though.

--- Reading in two triangulations

     std::ifstream in;

     in.open("input/solid.inp");
  GridIn<3> grid_in_sld;
  grid_in_sld.attach_triangulation(triangulation_sld);
  grid_in_sld.read_abaqus(in);

  std::ifstream in2;
  in2.open("input/flat3x3.inp");
  GridIn<2,3> grid_in_sh;
  grid_in_sh.attach_triangulation(triangulation_sh);
  grid_in_sh.read_abaqus(in2);

--- Partitioning the two meshes

    GridTools::partition_triangulation_zorder(
          Utilities::MPI::n_mpi_processes(mpi_communicator), triangulation_sld); // maybe divide n_mpi by 2
    GridTools::partition_multigrid_levels(triangulation_sld);
    auto construction_data_sld =
            TriangulationDescription::Utilities::create_description_from_triangulation(
                  triangulation_sld,  mpi_communicator,
                  TriangulationDescription::Settings::construct_multigrid_hierarchy);
    tria_pft_sld.create_triangulation(construction_data_sld);

    GridTools::partition_triangulation_zorder(
          Utilities::MPI::n_mpi_processes(mpi_communicator), triangulation_sh);// maybe divide n_mpi by 2
    GridTools::partition_multigrid_levels(triangulation_sh);
    auto construction_data_sh =
            TriangulationDescription::Utilities::create_description_from_triangulation(
                  triangulation_sh,  mpi_communicator,
                  TriangulationDescription::Settings::construct_multigrid_hierarchy);
    tria_pft_sh.create_triangulation(construction_data_sh);

If I do this, then it seems that I'd end up with two conflicting partitions.  If I have n_mpi_processes=8 , for example, this approach would assign ~1/8th of the solid mesh to process 1 AND ~1/8th of the shell mesh to process 1.  This seems like a problem.

I was thinking that maybe I could set it up such that the larger solid mesh could be split into 6 partitions over processes 0-5.  Then the smaller shell mesh could be split into 2 partitions using processes 6 and 7.  However, I didn't see any way to have the process numbers start at a non-zero number.

As you pointed out, I will have an issue with the off-diagonal matrices and the coupling between the two triangulations.  Is it possible for me to build up all of the sparsity patterns and matrices and then perform the partitioning?

Thanks,

Alex

[Wolfgang Bangerth]

On 1/17/24 14:40, Alex Quinlan wrote:
>
> If I do this, then it seems that I'd end up with two conflicting
> partitions.  If I have n_mpi_processes=8 , for example, this approach
> would assign ~1/8th of the solid mesh to process 1 AND ~1/8th of the
> shell mesh to process 1.  This seems like a problem.

No, this is exactly what you want. That's because you will have phases
in your program when you work on the solid mesh, and you want all 8
processes to participate in that. And then you have a phase where you
going to do something on the surface mesh and, again, you want to have
all 8 processes work on that.

It is *so* much harder to write software in which different processes do
different things at the same time. Don't go there.

> As you pointed out, I will have an issue with the off-diagonal matrices
> and the coupling between the two triangulations.  Is it possible for me

> to build up all of the sparsity patterns and matrices /and then /perform
> the partitioning?

You could of course, but then what's the point of using multiple
processors if you do the bulk of the work on only one process?

Best
W.

[Alex Quinlan]

Thanks for you response, Wolfgang.  The triangulation partitioning makes a lot of sense now.

In that case, perhaps my problem is not too difficult.  I am getting caught up with the PETScWrappers::MPI::BlockSparseMatrix.reinit()  method.

This is my approach:

Distributing dofs and collecting the locally owned and locally relevant sets.

// Solid

dof_handler_sld.distribute_dofs(fe_sld);

sld_owned_dofs = dof_handler_sld.locally_owned_dofs();

sld_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler_sld);

// Shell

dof_handler_sh.distribute_dofs(fe_sh);

sh_owned_dofs = dof_handler_sh.locally_owned_dofs();

sh_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler_sh);

Then I apply boundary conditions to both solid and shell dof_handlers. Then I define the coupling points and retrieve the coupled DOFs:

sld_coup_dofs = nodal_coupling(coup_points , dof_handler_sld, tria_pft_sld.n_vertices() );

sh_coup_dofs = nodal_coupling(coup_points , dof_handler_sh, tria_pft_sh.n_vertices() );

const std::vector<unsigned int> dofs_per_block = {n_dof_sld, n_dof_sh};

const std::vector<unsigned int> locally_owned_sizes = {sld_owned_dofs.size(), sh_owned_dofs.size()};

I reinitialize each of the blocks in my sparsity pattern:

BlockDynamicSparsityPattern dsp(dofs_per_block , dofs_per_block);

for (unsigned int i=0; i<dofs_per_block.size(); ++i)

  for (unsigned int j=0; j<dofs_per_block.size(); ++j)

      dsp.block(i,j).reinit(dofs_per_block[i], dofs_per_block[j]);

Then I use the solid and shell dof handlers to set up sub-matrices A and D.

dsp.collect_sizes();

DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); // A

DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1)); // D

Then, for the off-diagonal sub-matrices B and D, I add the coupled DOFs to the sparsity pattern.

for (unsigned int s = 0; s < sld_coup_dofs.size() ; s++)

  for (unsigned int d = 0; d < 3 ; d++)   {

    dsp.block(0,1).add(sld_coup_dofs[s][d] , sh_coup_dofs[s][d] ); // B

    dsp.block(1,0).add(sh_coup_dofs[s][d] , sld_coup_dofs[s][d] ); // C

  }

Where I start to get confused is when I try to reinitialize the system matrix.  What I want to do is something like:

system_matrix.reinit(dsp, mpi_communicator);

This has my MPI communicator and the sparsity pattern that I’ve built up. However, this isn’t a valid call to PETScWrappers::MPI::BlockSparseMatrix.reinit(). There’s a similar function that takes the arguments ( const std::vector< IndexSet > &sizes, const BlockDynamicSparsityPattern &bdsp, const MPI_Comm com )

I don’t really understand what I would put in for the "sizes" vector. What exactly am I trying to pass with this argument? Is it all of the locally owned/relevant dofs? Do I just combine the vector of locally owned shell dofs and locally owned solid dofs?

Thanks,

Alex

[Wolfgang Bangerth]

On 1/18/24 12:56, Alex Quinlan wrote:
> Where I start to get confused is when I try to reinitialize the system
> matrix.  What I want to do is something like:
>
> system_matrix.reinit(dsp, mpi_communicator);
>
>
> This has my MPI communicator and the sparsity pattern that I’ve built
> up. However, this isn’t a valid call to
> PETScWrappers::MPI::BlockSparseMatrix.reinit(). There’s a similar

> function that takes the arguments ( *const std::vector< IndexSet
> <https://www.dealii.org/current/doxygen/deal.II/classIndexSet.html> > &sizes,*const BlockDynamicSparsityPattern <https://www.dealii.org/current/doxygen/deal.II/classBlockDynamicSparsityPattern.html> &bdsp, const MPI_Comm <https://www.dealii.org/current/doxygen/deal.II/classMPI__Comm.html> com )

>
>
> I don’t really understand what I would put in for the "sizes" vector.
> What exactly am I trying to pass with this argument? Is it all of the
> locally owned/relevant dofs? Do I just combine the vector of locally
> owned shell dofs and locally owned solid dofs?

The argument is poorly named. It is a vector of index sets (in your case
a vector of 2 index sets) each of which contains which of the elements
of a block are locally owned. Perhaps the corresponding function in
other block classes uses a better named and/or documented argument. (The
name dates back to a time when we just passed down how many rows each
process owns, rather than being explicit -- with an index set -- which
rows these actually are.)

This would probably be worth fixing. Would you be willing to write a
patch that renames the argument?

Best
W.

[Alex Quinlan]

Thanks, Wolfgang.

I seem to have it working now, though I've hit some other unrelated snags that I need to resolve.  Once I get those fixed and confirm that my program is running correctly, I would be willing to work on a patch. 

Thanks,

Alex

[Wolfgang Bangerth]

On 1/24/24 09:07, Alex Quinlan wrote:
>
> I seem to have it working now, though I've hit some other unrelated
> snags that I need to resolve.  Once I get those fixed and confirm that
> my program is running correctly, I would be willing to work on a patch.

Excellent, thank you in advance!
Best
W.

[Alex Quinlan]

I am having more issues with MPI and the Block Matrix.  This time, it is the line:

system_matrix.compress(VectorOperation::add);

I get a PETSc error when I do this.  The error occurs within petsc_matrix_base.cc:275

// flush buffers
PetscErrorCode ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
AssertThrow(ierr == 0, ExcPETScError(ierr));

I've been going through this with GDB and found that one processor will have no problem with line, while all of the others create an error that looks like this:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at global row/column (36, 48) into matrix
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.12.4, Feb, 04, 2020

I wanted to see if there was anything special about the row/column (36,48), so I've overlapped an image of the sparsity patterns for each CPU.

The red image came from:

if (this_mpi_process == 0)    {
      sparsity_pattern.copy_from(dsp);
      sparsity_pattern.print_svg(out4);
    }

while the blue image came from this_mpi_process == 1.   I've marked the pixel with this matrix-entry in black in an image. It's a little hard to see, but it's just first entry into the blue, on the border with the purple overlapped area.  Therefore, it is not included in the sparsity pattern for cpu0.  I expect to get an error if cpu0 tries to do anything with this entry, but I don't understand why the system_matrix.compress(VectorOperation::add) is triggering this?  I wasn't able to find an answer in the glossary entry https://www.dealii.org/current/doxygen/deal.II/DEALGlossary.html#GlossCompress .

So, it seems that I'm doing something wrong with either my sparsity pattern or system matrix or maybe the constraints such that I can't properly do the system_matrix.compress() operation.  I had this working in serial, but this problem arose when I switched to MPI.

My setup and assembly looks like this:

void ElasticProblem::setup_system()
  {

    // --- Create p:f:t's

    // Solid

    GridTools::partition_triangulation_zorder( Utilities::MPI::n_mpi_processes(mpi_communicator), triangulation_sld); // maybe divide n_mpi by 2

    GridTools::partition_multigrid_levels(triangulation_sld);
    auto construction_data_sld = TriangulationDescription::Utilities::create_description_from_triangulation(
                  triangulation_sld,  mpi_communicator,  TriangulationDescription::Settings::construct_multigrid_hierarchy);

    tria_pft_sld.create_triangulation(construction_data_sld);

    // Shell
    GridTools::partition_triangulation_zorder(Utilities::MPI::n_mpi_processes(mpi_communicator), triangulation_sh);

    GridTools::partition_multigrid_levels(triangulation_sh);
    auto construction_data_sh = TriangulationDescription::Utilities::create_description_from_triangulation(
                  triangulation_sh, mpi_communicator, TriangulationDescription::Settings::construct_multigrid_hierarchy);
    tria_pft_sh.create_triangulation(construction_data_sh);

    // Solid
    dof_handler_sld.distribute_dofs(fe_sld);
    sld_owned_dofs = dof_handler_sld.locally_owned_dofs();
    sld_relevant_dofs =  DoFTools::extract_locally_relevant_dofs(dof_handler_sld);

    // Shell
    dof_handler_sh.distribute_dofs(fe_sh);
    sh_owned_dofs = dof_handler_sh.locally_owned_dofs();
    sh_relevant_dofs =    DoFTools::extract_locally_relevant_dofs(dof_handler_sh);

    const unsigned int n_dof_sld = dof_handler_sld.n_dofs();

    const unsigned int n_dof_sh  = dof_handler_sh.n_dofs();

    std::string solidbcfile = "input/solid.csv";  
    std::string shellbcfile = "input/shell.csv";

    constraints_sld.clear();
    nodal_bcs<3,3>(solidbcfile,  dof_handler_sld   ,constraints_sld);
    constraints_sld.close();
   

    constraints_sh.clear();
    nodal_bcs<2,3>(shellbcfile,  dof_handler_sh   ,constraints_sh);
    constraints_sh.close();

    constraints.copy_from(constraints_sld);
    constraints_sh.shift(n_dof_sld);
    constraints.merge(constraints_sh);  

    const std::vector<unsigned int> dofs_per_block = {n_dof_sld, n_dof_sh};
    const std::vector<unsigned int> locally_owned_sizes = {sld_owned_dofs.size(), sh_owned_dofs.size()};

    const std::vector<IndexSet> v_locown = {sld_owned_dofs, sh_owned_dofs};
    const std::vector<IndexSet> v_locrelv = {sld_relevant_dofs, sh_relevant_dofs};


    BlockDynamicSparsityPattern dsp(v_locrelv);

    dsp.collect_sizes();
    DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); // A

    DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1)); // D

    SparsityTools::distribute_sparsity_pattern(dsp, v_locown, mpi_communicator, v_locrelv);

    system_matrix.reinit(v_locown, dsp, mpi_communicator);
    system_rhs.reinit(v_locown, mpi_communicator);
    solution.reinit(dofs_per_block);
    system_matrix.collect_sizes();
}

void ElasticProblem::assemble_system()
  {
system_rhs    = 0;
system_matrix = 0;

... Create FE values objects ...

for (const auto &cell : dof_handler_sld.active_cell_iterators())

          if (cell->is_locally_owned())
            {
        cell_matrix_sld = 0;
         cell_rhs_sld    = 0;

 for (const unsigned int i : fe_values_sld.dof_indices())
       for (const unsigned int j : fe_values_sld.dof_indices())
         for (const unsigned int q_point : fe_values_sld.quadrature_point_indices())
        {
            cell_matrix_sld(i,j) = ...long equation...
        }
     
 cell->get_dof_indices(local_dof_indices_sld);

 constraints.distribute_local_to_global( cell_matrix_sld, cell_rhs_sld, 

                          local_dof_indices_sld, system_matrix.block(0,0), system_rhs.block(0));

    }  // end of cell loop

MPI_Barrier(mpi_communicator);
system_matrix.compress(VectorOperation::add);  // <---- This is where my issue is
system_rhs.compress(VectorOperation::add);

  // Onto Shell element assembly ...
...
...
}

The same error occurs if I skip the system_matrix.compress() at the end of the solid elements, and instead wait until the end of the shell element assembly.

I've looked at Step-55 and Step-70, which use the MPI::BlockSparseMatrix, but I don't see what I'm missing.  These examples use an index set locally_relevant_dofs_per_processor, but I am using a vector of index sets, since I have two DoF handlers. 

Can anyone see where I may have gone wrong or what I have omitted? 

Many thanks,

Alex

[Wolfgang Bangerth]

On 2/1/24 14:50, Alex Quinlan wrote:
>
> Can anyone see where I may have gone wrong or what I have omitted?

Alex:
I cannot see what the issue is, and I doubt anyone can without spending a
substantial amount of time debugging. I am slightly confused, though, by the
error message. I *believe* that what it is saying is that you are trying to
write into a matrix entry that is not part of the matrix -- you can only write
into matrix entries you have previously said exist when you passed the
sparsity pattern to the matrix. So what you need to find out is who owns row
36, and whether the sparsity pattern you provided on that process included
entry (36,48). My best guess is that that entry was not present in the
sparsity pattern on the process that owns row 36, and then your goal needs to
be to figure out why not -- in other words, why you did not foresee that the
entry that you ultimately end up writing into was not included when you
constructed the sparsity pattern.

[Alex Quinlan]

Thanks, Wolfgang.  Sorry for the big code dump.  I've been trying to distill the problem down to find the issue.  I've now eliminated the shell elements and the solid-shell coupling. In doing so, I've found a strange change to my sparsity pattern.  The overlap of Process0 and Process1 is not symmetric. The PETSc error occurs in this area where I'd expect to see symmetry (at entry 36,48).

The geometry of the problem is shown in the picture below.  I've colored the cells and vertices based on their ownership.

The dynamic sparsity pattern comes from this approach:

    dof_handler_sld.distribute_dofs(fe_sld);
    sld_owned_dofs = dof_handler_sld.locally_owned_dofs();

    DynamicSparsityPattern debug_dsp(sld_owned_dofs);
    DoFTools::make_sparsity_pattern(dof_handler_sld, debug_dsp);

It makes sense that I should have an entry at position (36,48), since dof 36 and dof 48 share a cell.  Therefore, I think my stiffness matrix process is okay.  My problem seems to be with my sparsity pattern.

So, I'm confused as to why Proc0 would not have any Sparsity Pattern entries for (36,48) when Proc1 does indeed have (48,36).

Any thoughts on where I might be causing this problem?

Best regards,

Alex

[Wolfgang Bangerth]

On 2/5/24 12:26, Alex Quinlan wrote:
>
> The dynamic sparsity pattern comes from this approach:
>
>     dof_handler_sld.distribute_dofs(fe_sld);
>     sld_owned_dofs = dof_handler_sld.locally_owned_dofs();
>
>     DynamicSparsityPattern debug_dsp(sld_owned_dofs);
>     DoFTools::make_sparsity_pattern(dof_handler_sld, debug_dsp);
>
> It makes sense that I should have an entry at position (36,48), since
> dof 36 and dof 48 share a cell.  Therefore, I think my stiffness matrix
> process is okay.  My problem seems to be with my sparsity pattern.
> So, I'm confused as to why Proc0 would not have any Sparsity Pattern
> entries for (36,48) when Proc1 does indeed have (48,36).

I don't know whether that's going to fix your problem, but the canonical
place to look up how to generate sparsity patterns in parallel is
step-40, where the corresponding code looks like this:

DynamicSparsityPattern dsp(locally_relevant_dofs);

DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
SparsityTools::distribute_sparsity_pattern(dsp,
dof_handler.locally_owned_dofs(),
mpi_communicator,
locally_relevant_dofs);

Best
W.