Periodic bc demos fail with Trilinos

[Garth Wells]

Error is thrown when run. Likely because of insertion into non-preallocated matrix position. Could also explain why periodic demos have been reported to be very slow.

[Anders Logg (logg)]

On Sun, Oct 02, 2011 at 05:02:48PM -0000, Launchpad Bug Tracker wrote:
> Garth Wells (garth-wells) has assigned this bug to you for DOLFIN:
>
> Error is thrown when run. Likely because of insertion into non-
> preallocated matrix position. Could also explain why periodic demos have
> been reported to be very slow.
>
> ** Affects: dolfin
> Importance: Medium
> Assignee: Anders Logg (logg)
> Status: Confirmed

Setting periodic bcs reduces the number of nonzeros on a row from n to
2 where n is likely > 2, so I would assume that would be ok.

Could it instead be modifying the matrix after apply has been called,
or (for PETSc) switching between add and insert without calling flush?

--
Anders

[Andre Massing (massing)]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I looked into this and it seems to me that is a bit more involved, so
here is what I have come up with so far. comments, corrections etc.
are very welcome!

At you might have seen I stumbled across 2 other bugs 891448 and 891452
which did not make it easier to find reasons or solutions for that bug :)

The implementation of PeriodicBC implicitly assumes that master and
slave dofs have the same row sparsity pattern as revealed by the
following lines ( PeriodicBC.cpp:180 )

 // Add slave rows to master rows
  for (uint i = 0; i < num_dof_pairs; ++i)
  {
    // Add slave row to master row in A
    if (A)
    {
      std::vector<uint> columns;
      std::vector<double> values;
      A->getrow(slave_dofs[i], columns, values);
      A->add(&values[0], 1, &master_dofs[i], columns.size(), &columns[0]);
      A->apply("add");
    }
However, it cannot be assumed that the row sparsity patter for the
master and slave are the same, since the sparsity pattern builder does
not take PeriodicBC into account (for instance by adding connectivity
information for dofs which are mapped to each other)

In addition states trilinos documentation that

"""
Note that, even after a matrix is constructed (FillComplete has been
called), it is possible to update existing matrix entries. It is not
possible to create new entries.
"""
see
http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/classEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5

And the program crashes exactly when a row patterns occur that do not
match (So a crossed mesh would be actually a "workaround" ;) ),
since the matrix is finalized when passed to the PeriodicBC::apply
function.

One could try to circumvent this by finalizing the matrix in the solve
function *after* assembling and the call to all boundary conditions,
which would require to
1. fixed the other aforementioned bugs regarding the unintentional
finalization of the Matrix
2. introduce a way to do several boundary condition apply calls
without assembling the matrix (up to now, apply does that in any case)

Nevertheless it still feels wrong, since the main reason for the bug
is the incomplete sparsity pattern which will result in bad
performance anyway. Or to put in a very "handwaving way", since the
Functionspace does not have a notion of constraints and we therefore
follow the procedure 1. finalize assemble 2. then modify A, b on a
algebraic level to take constraints into account, constraints are not
considered when the sparsity has to be computed.

Sorry for the detailed and elaborate email, I was to tired to write a
short one :)

Andre

On 11/16/2011 04:45 PM, Andre Massing wrote:
> ** Changed in: dolfin Status: Confirmed => In Progress
>
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.17 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJOxHH4AAoJEA79ggnbq9dmS+8H/jwNcIukNd8lgRcTE+WkNiGR
/iz3uEPSEolRf7AOdf/jOhl6jqz1rpwuRNI2ZR7wGPJbC3vM3fiNQClBw9A5pKOv
shTOHKxzLpZx/wv3pJz4lcNScBcSL0VdWbzTRMTOg/QY8LXbr5MpsGniZnOz1bGw
ZzUVyn1ZU4/N9oV7zX1qjOrN+m0AEwXxCMk3120CftqidaxSMYrTsuWuM/ndAWaa
EaM3dXi1sem1hF5OeJ3Bd7jNvaeZyS5dhcQ...

[Garth Wells (garth-wells)]

On 17 Nov 2011, at 02:31, Andre Massing <email address hidden> wrote:

> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> I looked into this and it seems to me that is a bit more involved, so
> here is what I have come up with so far. comments, corrections etc.
> are very welcome!
>
> At you might have seen I stumbled across 2 other bugs 891448 and 891452
> which did not make it easier to find reasons or solutions for that bug :)
>
> The implementation of PeriodicBC implicitly assumes that master and
> slave dofs have the same row sparsity pattern as revealed by the
> following lines ( PeriodicBC.cpp:180 )
>
> // Add slave rows to master rows
> for (uint i = 0; i < num_dof_pairs; ++i)
> {
> // Add slave row to master row in A
> if (A)
> {
> std::vector<uint> columns;
> std::vector<double> values;
> A->getrow(slave_dofs[i], columns, values);
> A->add(&values[0], 1, &master_dofs[i], columns.size(), &columns[0]);
> A->apply("add");
> }
> However, it cannot be assumed that the row sparsity patter for the
> master and slave are the same, since the sparsity pattern builder does
> not take PeriodicBC into account (for instance by adding connectivity
> information for dofs which are mapped to each other)
>
> In addition states trilinos documentation that
>
> """
> Note that, even after a matrix is constructed (FillComplete has been
> called), it is possible to update existing matrix entries. It is not
> possible to create new entries.
> """
> see
> http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/classEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
>
> And the program crashes exactly when a row patterns occur that do not
> match (So a crossed mesh would be actually a "workaround" ;) ),
> since the matrix is finalized when passed to the PeriodicBC::apply
> function.
>
> One could try to circumvent this by finalizing the matrix in the solve
> function *after* assembling and the call to all boundary conditions,
> which would require to
> 1. fixed the other aforementioned bugs regarding the unintentional
> finalization of the Matrix
> 2. introduce a way to do several boundary condition apply calls
> without assembling the matrix (up to now, apply does that in any case)
>
> Nevertheless it still feels wrong, since the main reason for the bug
> is the incomplete sparsity pattern which will result in bad
> performance anyway. Or to put in a very "handwaving way", since the
> Functionspace does not have a notion of constraints and we therefore
> follow the procedure 1. finalize assemble 2. then modify A, b on a
> algebraic level to take constraints into account, constraints are not
> considered when the sparsity has to be computed.
>

I registered a blueprint a while ago on this, which proposed taking the bcs into account when building the sparsity pattern. Modifying the sparsity of a matrix post-assembly is bad. Some backends allow it, but performance is usually poor, and other backends don't permit it.

Garth

> Sorry for the detailed and elaborate email, I was to tired to write a
> short one :)
>
> Andre
>
>
>
>
>
>
>
> On 11/16/2011 04:45 PM, Andre Massing wrot...

[Anders Logg (logg)]

On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> I looked into this and it seems to me that is a bit more involved, so
> here is what I have come up with so far. comments, corrections etc.
> are very welcome!
>
> At you might have seen I stumbled across 2 other bugs 891448 and 891452
> which did not make it easier to find reasons or solutions for that bug :)
>
> The implementation of PeriodicBC implicitly assumes that master and
> slave dofs have the same row sparsity pattern as revealed by the
> following lines ( PeriodicBC.cpp:180 )
>
> // Add slave rows to master rows
> for (uint i = 0; i < num_dof_pairs; ++i)
> {
> // Add slave row to master row in A
> if (A)
> {
> std::vector<uint> columns;
> std::vector<double> values;
> A->getrow(slave_dofs[i], columns, values);
> A->add(&values[0], 1, &master_dofs[i], columns.size(), &columns[0]);
> A->apply("add");
> }
> However, it cannot be assumed that the row sparsity patter for the
> master and slave are the same, since the sparsity pattern builder does
> not take PeriodicBC into account (for instance by adding connectivity
> information for dofs which are mapped to each other)
>
> In addition states trilinos documentation that
>
> """
> Note that, even after a matrix is constructed (FillComplete has been
> called), it is possible to update existing matrix entries. It is not
> possible to create new entries.
> """
> see
> http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/classEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
>
> And the program crashes exactly when a row patterns occur that do not
> match (So a crossed mesh would be actually a "workaround" ;) ),
> since the matrix is finalized when passed to the PeriodicBC::apply
> function.
>
> One could try to circumvent this by finalizing the matrix in the solve
> function *after* assembling and the call to all boundary conditions,
> which would require to
> 1. fixed the other aforementioned bugs regarding the unintentional
> finalization of the Matrix
> 2. introduce a way to do several boundary condition apply calls
> without assembling the matrix (up to now, apply does that in any case)
>
> Nevertheless it still feels wrong, since the main reason for the bug
> is the incomplete sparsity pattern which will result in bad
> performance anyway. Or to put in a very "handwaving way", since the
> Functionspace does not have a notion of constraints and we therefore
> follow the procedure 1. finalize assemble 2. then modify A, b on a
> algebraic level to take constraints into account, constraints are not
> considered when the sparsity has to be computed.
>
> Sorry for the detailed and elaborate email, I was to tired to write a
> short one :)

This may call for the introduction of "constrained function spaces",
being able to do things like

  V.constrain(bc0);
  V.constrain(bc1);

Another thing we have discussed is function spaces restricted to subdomains:

  V.restrict(subdomain)

I think this is an important feature to add for 1.1.

--
Anders

>
>
>
>
>
>
> On 11/16/2011 04:45 PM, Andre Massing wrote...

[Anders Logg (logg)]

Blueprint has been registered and this bug has been retargeted to trunk:

https://blueprints.launchpad.net/dolfin/+spec/constrained-functionspaces

[Andre Massing (massing)]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 11/17/2011 01:45 PM, Anders Logg wrote:
> On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote: I
> looked into this and it seems to me that is a bit more involved,
> so here is what I have come up with so far. comments, corrections
> etc. are very welcome!
>
> At you might have seen I stumbled across 2 other bugs 891448 and
> 891452 which did not make it easier to find reasons or solutions
> for that bug :)
>
> The implementation of PeriodicBC implicitly assumes that master
> and slave dofs have the same row sparsity pattern as revealed by
> the following lines ( PeriodicBC.cpp:180 )
>
> // Add slave rows to master rows for (uint i = 0; i <
> num_dof_pairs; ++i) { // Add slave row to master row in A if (A) {
> std::vector<uint> columns; std::vector<double> values;
> A->getrow(slave_dofs[i], columns, values); A->add(&values[0], 1,
> &master_dofs[i], columns.size(), &columns[0]); A->apply("add"); }
> However, it cannot be assumed that the row sparsity patter for the
> master and slave are the same, since the sparsity pattern builder
> does not take PeriodicBC into account (for instance by adding
> connectivity information for dofs which are mapped to each other)
>
> In addition states trilinos documentation that
>
> """ Note that, even after a matrix is constructed (FillComplete has
> been called), it is possible to update existing matrix entries. It
> is not possible to create new entries. """ see
> http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/classEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
>
> And the program crashes exactly when a row patterns occur that do
> not match (So a crossed mesh would be actually a "workaround" ;)
> ), since the matrix is finalized when passed to the
> PeriodicBC::apply function.
>
> One could try to circumvent this by finalizing the matrix in the
> solve function *after* assembling and the call to all boundary
> conditions, which would require to 1. fixed the other
> aforementioned bugs regarding the unintentional finalization of the
> Matrix 2. introduce a way to do several boundary condition apply
> calls without assembling the matrix (up to now, apply does that in
> any case)
>
> Nevertheless it still feels wrong, since the main reason for the
> bug is the incomplete sparsity pattern which will result in bad
> performance anyway. Or to put in a very "handwaving way", since
> the Functionspace does not have a notion of constraints and we
> therefore follow the procedure 1. finalize assemble 2. then modify
> A, b on a algebraic level to take constraints into account,
> constraints are not considered when the sparsity has to be
> computed.
>
> Sorry for the detailed and elaborate email, I was to tired to write
> a short one :)
>
>> This may call for the introduction of "constrained function
>> spaces", being able to do things like
>
>> V.constrain(bc0); V.constrain(bc1);
>
>> Another thing we have discussed is function spaces restricted to
>> subdomains:
>
>> V.restrict(subdomain)
>
>> I think this is an important feature to add for 1.1.

Yes, agree, some general notation of constraints would...

[Anders Logg (logg)]

On Thu, Nov 17, 2011 at 04:15:04PM +0100, Andre Massing wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
>
>
> On 11/17/2011 01:45 PM, Anders Logg wrote:
> > On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote: I
> > looked into this and it seems to me that is a bit more involved,
> > so here is what I have come up with so far. comments, corrections
> > etc. are very welcome!
> >
> > At you might have seen I stumbled across 2 other bugs 891448 and
> > 891452 which did not make it easier to find reasons or solutions
> > for that bug :)
> >
> > The implementation of PeriodicBC implicitly assumes that master
> > and slave dofs have the same row sparsity pattern as revealed by
> > the following lines ( PeriodicBC.cpp:180 )
> >
> > // Add slave rows to master rows for (uint i = 0; i <
> > num_dof_pairs; ++i) { // Add slave row to master row in A if (A) {
> > std::vector<uint> columns; std::vector<double> values;
> > A->getrow(slave_dofs[i], columns, values); A->add(&values[0], 1,
> > &master_dofs[i], columns.size(), &columns[0]); A->apply("add"); }
> > However, it cannot be assumed that the row sparsity patter for the
> > master and slave are the same, since the sparsity pattern builder
> > does not take PeriodicBC into account (for instance by adding
> > connectivity information for dofs which are mapped to each other)
> >
> > In addition states trilinos documentation that
> >
> > """ Note that, even after a matrix is constructed (FillComplete has
> > been called), it is possible to update existing matrix entries. It
> > is not possible to create new entries. """ see
> > http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/classEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
> >
> > And the program crashes exactly when a row patterns occur that do
> > not match (So a crossed mesh would be actually a "workaround" ;)
> > ), since the matrix is finalized when passed to the
> > PeriodicBC::apply function.
> >
> > One could try to circumvent this by finalizing the matrix in the
> > solve function *after* assembling and the call to all boundary
> > conditions, which would require to 1. fixed the other
> > aforementioned bugs regarding the unintentional finalization of the
> > Matrix 2. introduce a way to do several boundary condition apply
> > calls without assembling the matrix (up to now, apply does that in
> > any case)
> >
> > Nevertheless it still feels wrong, since the main reason for the
> > bug is the incomplete sparsity pattern which will result in bad
> > performance anyway. Or to put in a very "handwaving way", since
> > the Functionspace does not have a notion of constraints and we
> > therefore follow the procedure 1. finalize assemble 2. then modify
> > A, b on a algebraic level to take constraints into account,
> > constraints are not considered when the sparsity has to be
> > computed.
> >
> > Sorry for the detailed and elaborate email, I was to tired to write
> > a short one :)
> >
> >> This may call for the introduction of "constrained function
> >> spaces", being able to do things like
> >
> >> V.constrain(bc0); V.constrain(bc1);
> >
> >> Another thing we have discussed is function sp...

[Johan Hake (johan-hake)]

On Thursday November 17 2011 07:15:04 Andre Massing wrote:
> On 11/17/2011 01:45 PM, Anders Logg wrote:
> > On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote: I
> > looked into this and it seems to me that is a bit more involved,
> > so here is what I have come up with so far. comments, corrections
> > etc. are very welcome!
> >
> > At you might have seen I stumbled across 2 other bugs 891448 and
> > 891452 which did not make it easier to find reasons or solutions
> > for that bug :)
> >
> > The implementation of PeriodicBC implicitly assumes that master
> > and slave dofs have the same row sparsity pattern as revealed by
> > the following lines ( PeriodicBC.cpp:180 )
> >
> > // Add slave rows to master rows for (uint i = 0; i <
> > num_dof_pairs; ++i) { // Add slave row to master row in A if (A) {
> > std::vector<uint> columns; std::vector<double> values;
> > A->getrow(slave_dofs[i], columns, values); A->add(&values[0], 1,
> > &master_dofs[i], columns.size(), &columns[0]); A->apply("add"); }
> > However, it cannot be assumed that the row sparsity patter for the
> > master and slave are the same, since the sparsity pattern builder
> > does not take PeriodicBC into account (for instance by adding
> > connectivity information for dofs which are mapped to each other)
> >
> > In addition states trilinos documentation that
> >
> > """ Note that, even after a matrix is constructed (FillComplete has
> > been called), it is possible to update existing matrix entries. It
> > is not possible to create new entries. """ see
> > http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/c
> > lassEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
> >
> > And the program crashes exactly when a row patterns occur that do
> >
> > not match (So a crossed mesh would be actually a "workaround" ;)
> > ), since the matrix is finalized when passed to the
> > PeriodicBC::apply function.
> >
> > One could try to circumvent this by finalizing the matrix in the
> > solve function *after* assembling and the call to all boundary
> > conditions, which would require to 1. fixed the other
> > aforementioned bugs regarding the unintentional finalization of the
> > Matrix 2. introduce a way to do several boundary condition apply
> > calls without assembling the matrix (up to now, apply does that in
> > any case)
> >
> > Nevertheless it still feels wrong, since the main reason for the
> > bug is the incomplete sparsity pattern which will result in bad
> > performance anyway. Or to put in a very "handwaving way", since
> > the Functionspace does not have a notion of constraints and we
> > therefore follow the procedure 1. finalize assemble 2. then modify
> > A, b on a algebraic level to take constraints into account,
> > constraints are not considered when the sparsity has to be
> > computed.
> >
> > Sorry for the detailed and elaborate email, I was to tired to write
> > a short one :)
> >
> >> This may call for the introduction of "constrained function
> >> spaces", being able to do things like
> >>
> >> V.constrain(bc0); V.constrain(bc1);
> >>
> >> Another thing we have discussed is function spaces restricted to
> >> subdomains:
...

[Andre Massing (massing)]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 11/17/2011 07:34 PM, Johan Hake wrote:
> On Thursday November 17 2011 07:15:04 Andre Massing wrote:
>> On 11/17/2011 01:45 PM, Anders Logg wrote:
>>> On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote:
>>> I looked into this and it seems to me that is a bit more
>>> involved, so here is what I have come up with so far. comments,
>>> corrections etc. are very welcome!
>>>
>>> At you might have seen I stumbled across 2 other bugs 891448
>>> and 891452 which did not make it easier to find reasons or
>>> solutions for that bug :)
>>>
>>> The implementation of PeriodicBC implicitly assumes that
>>> master and slave dofs have the same row sparsity pattern as
>>> revealed by the following lines ( PeriodicBC.cpp:180 )
>>>
>>> // Add slave rows to master rows for (uint i = 0; i <
>>> num_dof_pairs; ++i) { // Add slave row to master row in A if
>>> (A) { std::vector<uint> columns; std::vector<double> values;
>>> A->getrow(slave_dofs[i], columns, values); A->add(&values[0],
>>> 1, &master_dofs[i], columns.size(), &columns[0]);
>>> A->apply("add"); } However, it cannot be assumed that the row
>>> sparsity patter for the master and slave are the same, since
>>> the sparsity pattern builder does not take PeriodicBC into
>>> account (for instance by adding connectivity information for
>>> dofs which are mapped to each other)
>>>
>>> In addition states trilinos documentation that
>>>
>>> """ Note that, even after a matrix is constructed (FillComplete
>>> has been called), it is possible to update existing matrix
>>> entries. It is not possible to create new entries. """ see
>>> http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/c
>>>
>>>
lassEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
>>>
>>> And the program crashes exactly when a row patterns occur that
>>> do
>>>
>>> not match (So a crossed mesh would be actually a "workaround"
>>> ;) ), since the matrix is finalized when passed to the
>>> PeriodicBC::apply function.
>>>
>>> One could try to circumvent this by finalizing the matrix in
>>> the solve function *after* assembling and the call to all
>>> boundary conditions, which would require to 1. fixed the other
>>> aforementioned bugs regarding the unintentional finalization of
>>> the Matrix 2. introduce a way to do several boundary condition
>>> apply calls without assembling the matrix (up to now, apply
>>> does that in any case)
>>>
>>> Nevertheless it still feels wrong, since the main reason for
>>> the bug is the incomplete sparsity pattern which will result in
>>> bad performance anyway. Or to put in a very "handwaving way",
>>> since the Functionspace does not have a notion of constraints
>>> and we therefore follow the procedure 1. finalize assemble 2.
>>> then modify A, b on a algebraic level to take constraints into
>>> account, constraints are not considered when the sparsity has
>>> to be computed.
>>>
>>> Sorry for the detailed and elaborate email, I was to tired to
>>> write a short one :)
>>>
>>>> This may call for the introduction of "constrained function
>>>> spaces", being able to do things like
>>>>
>>>> V.constrain(bc0); ...

[Anders Logg (logg)]

On Thu, Nov 17, 2011 at 08:42:20PM -0000, Andre Massing wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
>
> On 11/17/2011 07:34 PM, Johan Hake wrote:
> > On Thursday November 17 2011 07:15:04 Andre Massing wrote:
> >> On 11/17/2011 01:45 PM, Anders Logg wrote:
> >>> On Thu, Nov 17, 2011 at 02:31:20AM -0000, Andre Massing wrote:
> >>> I looked into this and it seems to me that is a bit more
> >>> involved, so here is what I have come up with so far. comments,
> >>> corrections etc. are very welcome!
> >>>
> >>> At you might have seen I stumbled across 2 other bugs 891448
> >>> and 891452 which did not make it easier to find reasons or
> >>> solutions for that bug :)
> >>>
> >>> The implementation of PeriodicBC implicitly assumes that
> >>> master and slave dofs have the same row sparsity pattern as
> >>> revealed by the following lines ( PeriodicBC.cpp:180 )
> >>>
> >>> // Add slave rows to master rows for (uint i = 0; i <
> >>> num_dof_pairs; ++i) { // Add slave row to master row in A if
> >>> (A) { std::vector<uint> columns; std::vector<double> values;
> >>> A->getrow(slave_dofs[i], columns, values); A->add(&values[0],
> >>> 1, &master_dofs[i], columns.size(), &columns[0]);
> >>> A->apply("add"); } However, it cannot be assumed that the row
> >>> sparsity patter for the master and slave are the same, since
> >>> the sparsity pattern builder does not take PeriodicBC into
> >>> account (for instance by adding connectivity information for
> >>> dofs which are mapped to each other)
> >>>
> >>> In addition states trilinos documentation that
> >>>
> >>> """ Note that, even after a matrix is constructed (FillComplete
> >>> has been called), it is possible to update existing matrix
> >>> entries. It is not possible to create new entries. """ see
> >>> http://trilinos.sandia.gov/packages/docs/r10.8/packages/epetra/doc/html/c
> >>>
> >>>
> lassEpetra__CrsMatrix.html#a1291bfb498b94ac4f1563b48d11dfae5
> >>>
> >>> And the program crashes exactly when a row patterns occur that
> >>> do
> >>>
> >>> not match (So a crossed mesh would be actually a "workaround"
> >>> ;) ), since the matrix is finalized when passed to the
> >>> PeriodicBC::apply function.
> >>>
> >>> One could try to circumvent this by finalizing the matrix in
> >>> the solve function *after* assembling and the call to all
> >>> boundary conditions, which would require to 1. fixed the other
> >>> aforementioned bugs regarding the unintentional finalization of
> >>> the Matrix 2. introduce a way to do several boundary condition
> >>> apply calls without assembling the matrix (up to now, apply
> >>> does that in any case)
> >>>
> >>> Nevertheless it still feels wrong, since the main reason for
> >>> the bug is the incomplete sparsity pattern which will result in
> >>> bad performance anyway. Or to put in a very "handwaving way",
> >>> since the Functionspace does not have a notion of constraints
> >>> and we therefore follow the procedure 1. finalize assemble 2.
> >>> then modify A, b on a algebraic level to take constraints into
> >>> account, constraints are not considered when the sparsity has
> >>> to be computed.
> >>>
> >>> Sorry for the detailed and elaborate email, I was to t...