Hey, the SparsityPattern class uses optimize_diagonal=true as a default. This always adds the diagonal entries. You can fix that by doing i.e. .reinit(n,n,0,false).
-- Timo Heister http://num.math.uni-goettingen.de/~heister On Fri, Sep 24, 2010 at 5:37 PM, fname lname <[email protected]> wrote: > Thanks all, sorry, I will try to explain myself a little better. > > I have several fully uncoupled systems of equations that correspond each one > to a finite element approximation for least squares formulation to solve a > differential equation. > > This uncoupled systems are what I want to solve including all of them in the > same matrix to solve, that gives me a block sparsity pattern, but one on > which the only blocks that could have non-zero elements would be the blocks > on the full matrix diagonal (because they are uncoupled). > > I am trying to define this particular sparsity pattern but when I do I get > off-diagonal blocks that are themselves diagonal matrices instead of "null" > or "empty" blocks without any non-zero elements, as I wished them to be. The > blocks in the diagonal are fine. > > This is to say I want a sparsity pattern like this: > > *** > *** > *** > *** > *** > *** > *** > *** > *** > *** > *** > *** > > And I'm getting a sparsity pattern like this: > > **** * * > *** * * * > *** * * * > * **** * > * *** * * > **** * * > * * **** > * * *** * > * **** * > * * * *** > * * * *** > * * **** > > The code actually works, but it is working with non-zero element blocks > outside the diagonal, so it is being to expensive in memory when I solve > more than a hundred sub-systems. > > I have tried moving the "compress" statement, defining a "null" sparsity > pattern as Timo suggested, by the off-diagonal matrices are still not empty. > > Thanks again, the code I have until now is down here. > > /******************* > /CODE > /******************* > > GridGenerator::hyper_cube (triangulation, -5, 5); > > triangulation.begin_active()->face(0)->set_boundary_indicator(0); > triangulation.begin_active()->face(1)->set_boundary_indicator(1); > triangulation.begin_active()->face(2)->set_boundary_indicator(2); > triangulation.begin_active()->face(3)->set_boundary_indicator(3); > triangulation.begin_active()->face(4)->set_boundary_indicator(4); > triangulation.begin_active()->face(5)->set_boundary_indicator(5); > > triangulation.refine_global (2); > > dof_handler.distribute_dofs (fe); > > CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); > > // I have also tried here > // CompressedSparsityPattern null_sparsity(dof_handler.n_dofs()); > > DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); > > sparsity_pattern.reinit(nd,nd); > > for (int i=0; i<nd; ++i) { > for (int j=0; j<nd; ++j) { > if (i==j) { > sparsity_pattern.block(i,j).copy_from(c_sparsity); > } else { > > sparsity_pattern.block(i,j).reinit(dof_handler.n_dofs(),dof_handler.n_dofs(),0); > > // I have also tried here > // sparsity_pattern.block(i,j).copy_from(null_sparsity); > // instead of the reinit without and making make_sparsity_pattern to this > null_sparsity block > > } > } > } > > sparsity_pattern.compress(); > > sparsity_pattern.collect_sizes (); > > system_matrix.reinit (sparsity_pattern); > > // .... and the rest of the code. > > On Fri, Sep 24, 2010 at 9:54 AM, Guido Kanschat <[email protected]> wrote: >> >> On 9/23/2010 7:16 PM, fname lname wrote: >>> >>> Thanks Timo, >>> >>> I've tried that already, it gets me the same matrix as if I do: >>> >>> >>> sparsity_pattern.block(i,j).reinit(dof_handler.n_dofs(),dof_handler.n_dofs(),0); >> >> Can't you call sparsity_pattern.block(i,j).compress() right there? >> >> G >> _______________________________________________ >> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii > > _______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
