[deal.II] Bug in sparse_matrix.templates.h

[Michał Wichrowski]

Dear all, 
Matrix multiplication requires resulting matrix to be initialized. Even if 
rebuild_sparsity_C = true .

template 
template 
void
SparseMatrix::mmult (SparseMatrix   &C,
 const SparseMatrix &B,
 const Vector&V,
 const bool   
rebuild_sparsity_C) const
{
  const bool use_vector = V.size() == n() ? true : false;
  Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
  Assert (cols != 0, ExcNotInitialized());
  Assert (B.cols != 0, ExcNotInitialized());
  Assert (C.cols != 0, ExcNotInitialized()); //HERE!!

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Bug in sparse_matrix.templates.h

[Michał Wichrowski]

Thanks Wolfgang. 
I think it would be helpful if there was something about it in 
documentation. 
Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] CUDA without nVidia

[Michał Wichrowski]

Dear all, 
Have anyone managed to compile and *run* CUDA code on non-nVidia device?
I've found HIPify: https://github.com/ROCm-Developer-Tools/HIPIFY 
that translates CUDA code to HIP code that can be then compiled and ran on 
both AMD and nVidia GPUs. 

The thing is, I'm looking for a new laptop. Since I'm planning switching to 
GPU I want to be able to compile, run and test my code locally. Choosing a 
hardware with nVidia chip narrows the choice, moreover AMD CPUs  seem to 
come with good enough built-in GPUs.

 Some years ago I've tried GPUOcelot and did not managed to make it work.

Best,
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/846802b0-3b0a-4933-92e3-3f1efe720ef5n%40googlegroups.com.

[deal.II] Re: CUDA without nVidia

[Michał Wichrowski]

Michal,
>
> Support for AMD through HIP is very high on my todo list but the compiler 
> is not that great right now, i.e. the compiler can (and will) generate 
> buggy code.
>
Ok, so even if Deal.II will be ported I can assume so for other libraries.  
Good to know that HIP does  require porting.

I've tried to hipify the code a few months ago but it didn't go too well. 
> My plan is to work on support for AMD after we branch the release so that 
> we don't have to support buggy compilers. I am pretty confident that we 
> will have support for HIP by the end of the Summer maybe earlier.
>
Great news! 

>  
>
Best,
>
> Bruno
>
> Thanks!
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c1855148-6990-4d68-aed0-03adc90ea838n%40googlegroups.com.

[deal.II] deal.II-based FSI matrix-free solver / looking for a post-doc

[Michał Wichrowski]

Dear all,
 I have developed a parallel matrix-free monolithic FSI solver in the ALE 
frame of reference. Exemplary results [1 <https://youtu.be/7TaC1dTfAQY>, 2 
<https://youtu.be/TRnlt_eC32Q>], more details in my PhD thesis [3] 
<http://www.ippt.pan.pl/repository/doktoraty/open/2021wichrowski_m_doktorat.pdf>
.

*Highlights*:

-> The proposed monolithic predictor-corrector time integration scheme 
requires solving the generalized Stokes problem with discontinuous variable 
coefficients. This has to be done one or two times per time step, depending 
on the variant of the scheme.

-> I have introduced a new multilevel preconditioner, for a 3D case the 
time required to solve a system consisting of 421 million unknowns using 96 
cores is 178 seconds.

-> As demonstrated by the movie[1], I am able to solve 3D problems at 
moderate Reynolds numbers (time step size of 0.01 did not result in any 
stability issues with Re=2k). 

->The preconditioner from the thesis can still be greatly improved. I have 
a proof of concept demonstrating the basic properties of a new idea of 
smoother. The results look very promising and I'm curious how it would 
work. Moreover, the algorithm is suitable for GPU that could result in 
further speedup.
 
[1] 3D "Turek benchmark" at Re = 2000 (30M DoF)
https://youtu.be/7TaC1dTfAQY
[2] Water hammer in elastic tube (12M DoF)
https://youtu.be/TRnlt_eC32Q

Those results were obtained as a part of my PhD thesis. I have finally 
managed to submit it in January, the defence will be held in 3 weeks 
(expect reviews next week).
I am looking for a post-doc (preferably deal.II-related). I'll be happy to 
get an opportunity to continue my work.

Best,
Michał Wichrowski

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/00edc1de-b7fd-4aa7-b0cc-709d66dabdf6n%40googlegroups.com.


re-pipe3D.mp4
Description: video/mp4

[deal.II] Re: deal.II-based FSI matrix-free solver / looking for a post-doc

[Michał Wichrowski]

>
> @Peter
Thanks!
@Bruno,
I've sent you an email last week.

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/153b6804-738b-44b6-8ef4-2d8e1b10655bn%40googlegroups.com.

[deal.II] Building p4est with no system-wide zlib

[Michał Wichrowski]

Dear all,
I am trying to set up deal.II on a computer with no system-wide zlib. I do 
not have sudo privileges, thus I cannot just install, so I compiled it from 
the source and installed in my home directory. Now I am trying to build 
p4est using it. According to the message thrown by the scrip 
(https://www.dealii.org/current/external-libs/p4est.html) I should 
export CPPFLAGS="-DSC_LOG_PRIORITY=SC_LP_ESSENTIAL  
-I/export/home/mwichro/lib/zlib/include/"
export LDFLAGS="-L/export/home/mwichro/lib/zlib/lib/"
Apparently setting those does not resolve the problem, the script still 
complains about zlib.

By going through what the script really does, I think the following line 
triggers an error:
grep -q 'P4EST_HAVE_ZLIB *1' "$BUILD_FAST/src/p4est_config.h" \ || bdie 
"$MISSING_ZLIB_MESSAGE"

As far as I understood, it looks into p4est-build/FAST/src/p4est_config.h . 
That file on my system does not exist. I also went through the config.log, 
it looks like it contains:
#define LIBS "  -lz -lm   "
so I guess it links with zlib?


By the way, I moved to Heidelberg and I am happy to continue my work with 
deal.II library.

Best,
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/31131e7f-0409-4889-9475-d2165fb26ae4n%40googlegroups.com.

Re: [deal.II] Building p4est with no system-wide zlib

[Michał Wichrowski]

@Wolfgang:
I commented out that lines, compiled deal.II and it looks OK, at least 
`make test` passes and  step-37 works in both release and debug mode. 

@Marc
I had the same idea, but there was nothing in `configure --help` to put 
custom zlib patch.

Thank you for the answers. For reference, I was using p4est-2.3.2 

piątek, 17 września 2021 o 08:36:09 UTC+2 mafe...@gmail.com napisał(a):

> Hi Michał,
>
> maybe it will help to invoke `configure --help` on the source of your 
> p4est tarball. Hopefully there will be instructions on how to connect 
> custom zlib libraries -- I couldn't find anything in their README :(
>
> Best,
> Marc
>
> On Thursday, September 16, 2021 at 6:26:23 PM UTC-6 Wolfgang Bangerth 
> wrote:
>
>>
>> > I am trying to set up deal.II on a computer with no system-wide zlib. I 
>> > do not have sudo privileges, thus I cannot just install, so I compiled 
>> > it from the source and installed in my home directory. Now I am trying 
>> > to build p4est using it. According to the message thrown by the scrip 
>> > (https://www.dealii.org/current/external-libs/p4est.html) I should 
>> > export CPPFLAGS="-DSC_LOG_PRIORITY=SC_LP_ESSENTIAL 
>> > -I/export/home/mwichro/lib/zlib/include/" 
>> > export LDFLAGS="-L/export/home/mwichro/lib/zlib/lib/" 
>> > Apparently setting those does not resolve the problem, the script still 
>> > complains about zlib. 
>> > 
>> > By going through what the script really does, I think the following 
>> line 
>> > triggers an error: 
>> > grep -q 'P4EST_HAVE_ZLIB *1' "$BUILD_FAST/src/p4est_config.h" \ || bdie 
>> > "$MISSING_ZLIB_MESSAGE" 
>> > 
>> > As far as I understood, it looks 
>> > into p4est-build/FAST/src/p4est_config.h . That file on my system does 
>> > not exist. I also went through the config.log, it looks like it 
>> contains: 
>> > #define LIBS "  -lz -lm   " 
>> > so I guess it links with zlib? 
>>
>> Michal: I don't know whether anyone still remembers how that 
>> installation script came about. What happens if you just remove the line 
>> that causes the script to abort? 
>>
>>
>> > By the way, I moved to Heidelberg and I am happy to continue my work 
>> > with deal.II library. 
>> Nice, glad to hear you landed in a good place a lot of us have good 
>> memories of! 
>>
>> Best 
>> Wolfgang 
>>
>> -- 
>>  
>> Wolfgang Bangerth email: 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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f49c5d4b-f9fa-40dc-ad12-c64481a32540n%40googlegroups.com.

[deal.II] Re: Elasticity Preconditioner

[Michał Wichrowski]

The simple answer is multigrid. But depending on how the problem looks like 
it will differ. If it is far from being incompressible standard MG will be 
fine. Trilinos AMG should also work, just take look at its settings (there 
is an option to specify that you are dealing with vector problem).
Best,
Michał


piątek, 24 września 2021 o 14:17:33 UTC+2 Konrad Simon napisał(a):

> Dear all, 
>
> A student and me are trying to deal with (parallel, static) linearized 
> elasticity similar to step-8.
> However, the problem is slightly different since the Lame parameters are 
> functions instead of constants. 
> Can anyone recommend good solvers & preconditioners for the resulting 
> system? Our boundary conditions fix all degrees of freedom in the kernel.
>
> Best,
> Konrad
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f9a5106b-cbba-40f5-8be7-b89c4473187en%40googlegroups.com.

[deal.II] PETSc autodetect

[Michał Wichrowski]

Dear all,
While I was setting up deal.II compilation on another machine I got stuck 
for weeks: it was crashing when running with MPI. It turned out that 
Deal.II managed to auto-detect PETSc (I did not turn it on, I am not using 
it), which was causing issues with the initialization of any MPI program. 

May I suggest changing Deal.II configuration script to disable 
autodetection of PETSc by default and do it only if the user explicitly 
asked for it? To be honest, I had a bad experience with it and I haven't 
found PETSc superior to Trilinos in any application.

Best, 
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0ade937d-3b3b-4816-a5a6-1ab7e00a9515n%40googlegroups.com.

[deal.II] PhD project opportunity at Heidelberg University

[Michał Wichrowski]

We (Guido Kanschat and I) offer a PhD project devising a code for 
fluid-structure interaction for cell movement, which should allow for 
experiments with various fluid and solid mechanics models in cooperation 
with experimentalists.

We are looking for somebody with experience with deal.II and preferably 
fluid or solid mechanics.

A master's degree is required when the person starts. Heidelberg University 
is an equal opportunity employer where women are particularly invited to 
apply.

If you are interested, don't hesitate to get in touch with me at 
wichrow...@uni-heidelberg.de

Best,
Michał Wichrowski

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2db1906d-a0a0-4899-95db-faef144468f8n%40googlegroups.com.

[deal.II] Block problem +MatrixFree

[Michał Wichrowski]

Dear all,
I'm dealing with Stokes problem that I want to solve using matrix-free 
methods. In tests I've alredy found  implementation of Stokes operator that 
uses block vectors. I'm goning to use multigrid method for saddle point 
problems similar to the one presented in:
An efficient smoother for the Stokes problem D Braess, R Sarazin - Applied 
Numerical Mathematics, 1997 - Elsevier

Let the Stokes matrix be:
(A B^T
B 0)

1) The method requires several multiplication by blocks A and B at each 
level. Is there any recommended way to implement it? ()

2) I've  found  MGTransferBlockMatrixFree class in developer version of 
deal.ii, however it assumes that each block have the same dof_hander (and 
thats for me is obviously not true). Is there existing MGTrasfrer for my 
problem?
BTW, why MGTrasfrer does work with several dof handres like MatrixFree?

3) Is there a way to use multigrid on one block? For example I would like 
to compute approximation of A^-1 using multigrid.
 

Michał Wichrowski

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Block problem +MatrixFree

[Michał Wichrowski]

W dniu wtorek, 26 września 2017 18:42:33 UTC+2 użytkownik Wolfgang Bangerth 
napisał:
>
> On 09/26/2017 09:42 AM, Michał Wichrowski wrote: 
> > 
> > 3) Is there a way to use multigrid on one block? For example I would 
> > like to compute approximation of A^-1 using multigrid. 
>
> This might be the way to go since applying multigrid on the whole 
> saddle-point problem is difficult (but feasible). 
>
> step-22 shows how to create a preconditioner that involves A^{-1} -- 
> take a particular look at the "Possibilities for extensions" section. I 
> think one of the Stokes or Navier-Stokes equation problems (54? 55? 56? 
> 57?) may even implement something that applies multigrid already to the 
> A block of the Stokes system. You may then want to replace that 
> multigrid method by a matrix-free version. I'll let others comment on 
> how this can best be done. 
>
> My approach is different, I'm using B-S method that I mentioned in my 
post. The smoother works on whole problem  (not on blocks), the smoothing 
operation is done by applying inverse of following matrix:
(S{A} B^T
B 0)
Where S(A) is some simple preconditioner for A (e.g. SSOR, Jaccobi etc).
I've already tested it for standard matrices and it works surprisingly well 
for standard Stokes equation. The only problem with doing it in matrix-free 
way is multiplication by block.

However in mine application slightly modified Stokes equation appears and I 
will need to do some multigid for blocks inside multigrid.
 

>
> Yours is a case that I think a lot of people would be interested in, 
> myself included. If you get it to work, could I interest you in thinking 
> about whether your program could be made into either a tutorial program 
> or a contribution to the code gallery? 
>

I'm probably  going to publish my method. Since I do not have mathematica 
proof for it I'm  thinking of  publishing source code. 

>
> Best 
>   W. 
>

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Block problem +MatrixFree

[Michał Wichrowski]

W dniu wtorek, 26 września 2017 20:37:59 UTC+2 użytkownik Daniel Arndt 
napisał:
>
> Michael,
>
> Adding to Wolfgang's response:
>
> Am Dienstag, 26. September 2017 17:42:41 UTC+2 schrieb Michał Wichrowski:
>>
>> Dear all,
>> I'm dealing with Stokes problem that I want to solve using matrix-free 
>> methods. In tests I've alredy found  implementation of Stokes operator that 
>> uses block vectors. I'm goning to use multigrid method for saddle point 
>> problems similar to the one presented in:
>> An efficient smoother for the Stokes problem D Braess, R Sarazin - 
>> Applied Numerical Mathematics, 1997 - Elsevier
>>
>> Let the Stokes matrix be:
>> (A B^T
>> B 0)
>>
>> 1) The method requires several multiplication by blocks A and B at each 
>> level. Is there any recommended way to implement it? ()
>>
> Have a look at tests/matrix_free/matrix_vector_stokes_base.cc for a 
> MatrixFree approach for a block system using BlockVectors.
> If you want the individual blocks A and B it might be best to write them 
> to implement them individually as derived classes from 
> MatrixFreeOperators::Base in a similar way as this 
> is done in the test mentioned above.
>
 

>  
>
>> 2) I've  found  MGTransferBlockMatrixFree class in developer version of 
>> deal.ii, however it assumes that each block have the same dof_hander (and 
>> thats for me is obviously not true). Is there existing MGTrasfrer for my 
>> problem?
>> BTW, why MGTrasfrer does work with several dof handres like MatrixFree?
>>
> No, this is currently not implemented but definitely on my TODO-list. This 
> would basically be combining MGTransferMatrixFree with the current 
> MGTransferBlockMatrixFree.
>  
>
I've started working on new version MGTransferBlockMatrixFree and I now see 
the problem. MGTransferBlockMatrixFree is easy to rewrite using vector of 
dof handres etc, but PreconditionMG requires MGTransfer that have 
 copy_to_mg and copy_from_mg that work on single dof handre. I am thinking 
of 3 way to resolve it:
1) Rewrite MGTransferBlockMatrixFree using  vector of dof handler, 
implement PreconditionBlockMG and probably Multigrid will also need block 
version
2) Use standard MGTransfer for individual blocks and implement 
PreconditionBlockMG and BlockMultigrid to work with several MGTransfer s.
3) Since MGTransferBlockMatrixFree requires dof handler in build (const 
DoFHandler< dim, dim > &mg_dof) method, the pointers to dof_handlers may be 
stored and used in copy_from_mg and copy_to_mg. 
I think 3) is simplest but also risky solution. 2) seems quite rational to 
me. 

>
> Best,
 Michał 

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Block problem +MatrixFree

[Michał Wichrowski]

>
>
> I've started working on new version MGTransferBlockMatrixFree and I now 
> see the problem. MGTransferBlockMatrixFree is easy to rewrite using vector 
> of dof handres etc, but PreconditionMG requires MGTransfer that have 
>  copy_to_mg and copy_from_mg that work on single dof handre. I am thinking 
> of 3 way to resolve it:
> 1) Rewrite MGTransferBlockMatrixFree using  vector of dof handler, 
> implement PreconditionBlockMG and probably Multigrid will also need block 
> version
> 2) Use standard MGTransfer for individual blocks and implement 
> PreconditionBlockMG and BlockMultigrid to work with several MGTransfer s.
> 3) Since MGTransferBlockMatrixFree requires dof handler in build (const 
> DoFHandler< dim, dim > &mg_dof) method, the pointers to dof_handlers may be 
> stored and used in copy_from_mg and copy_to_mg. 
> I think 3) is simplest but also risky solution. 2) seems quite rational to 
> me. 
>
>>
>> I've decided to implement approach 3). You can find it in attachment. 
Second file contains test on block Laplace problem. Main difference between 
this file and one of deal.II test are this lines:

  std::vector< const MGConstrainedDoFs * > mg_constrained_dofs_vector(nb, 
&mg_constrained_dofs );
  std::vector< const DoFHandler * > mg_dofs_vector(nb, &dof );
  MWichrowski::MGTransferBlockMatrixFree 
mg_transfer(mg_constrained_dofs_vector);
  mg_transfer.build(mg_dofs_vector);
 

> Best,
>  Michał 
>

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/*
 * MGTransferBlockMatrixFree.h
 *
 *  Created on: Sep 26, 2017
 *  Author: mwichro
 */

#ifndef MGTRANSFERBLOCKMATRIXFREE_H_
#define MGTRANSFERBLOCKMATRIXFREE_H_

#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace dealii;
namespace MWichrowski {
template 
class MGTransferBlockMatrixFree : public MGTransferBase>
{
public:
  /**
   * Constructor without constraint matrices. Use this constructor only with
   * discontinuous finite elements or with no local refinement.
   */
  MGTransferBlockMatrixFree (const unsigned int &n_b) //= default;
	:
	matrix_free_transfer(n_b),
	n_blocks(n_b)
{}

  /**
   * Constructor with constraints. Equivalent to the default constructor
   * followed by initialize_constraints().
   */
  MGTransferBlockMatrixFree (const std::vector< const MGConstrainedDoFs * > &mg_constrained_dofs);

  /**
   * Destructor.
   */
  virtual ~MGTransferBlockMatrixFree () = default;

  /**
   * Initialize the constraints to be used in build().
   */
  void initialize_constraints (const std::vector< const MGConstrainedDoFs* >  &mg_constrained_dofs);

  /**
   * Reset the object to the state it had right after the default constructor.
   */
  void clear ();


  /**
   * Actually build the information for the prolongation for each level.
   */
  void build (const std::vector * > &mg_dofs);

  /**
   * Prolongate a vector from level to_level-1 to level
   * to_level using the embedding matrices of the underlying finite
   * element. The previous content of dst is overwritten.
   *
   * @param to_level The index of the level to prolongate to, which is the
   * level of @p dst.
   *
   * @param src is a vector with as many elements as there are degrees of
   * freedom on the coarser level involved.
   *
   * @param dst has as many elements as there are degrees of freedom on the
   * finer level.
   */
  virtual void prolongate (const unsigned intto_level,
   LinearAlgebra::distributed::BlockVector   &dst,
   const LinearAlgebra::distributed::BlockVector &src) const;

  /**
   * Restrict a vector from level from_level to level
   * from_level-1 using the transpose operation of the prolongate()
   * method. If the region covered by cells on level from_level is
   * smaller than that of level from_level-1 (local refinement), then
   * some degrees of freedom in dst are active and will not be
   * altered. For the other degrees of freedom, the result of the restriction
   * is added.
   *
   * @param from_level The index of the level to restrict from, which is the
   * level of @p src.
   *
   * @param src is a vector with as many elements as there are degrees of
   * freedom on the finer level involved.
   *
   * @param dst has as many elements as there are degrees of freedom on the
   * coarser level.
   */
  virtual void restrict_and_add (const unsigned int from_level,
 LinearAlgebra::distributed::BlockVector   &dst,
 const LinearAlgebra::distributed::BlockV

[deal.II] Selective use of blocks in MatrixFree

[Michał Wichrowski]

Dear all,
I'm testing MatrixFree framework to build Stokes solver. I'm using Stokes 
operator copied from deal.ii tests.
1) I've found out that using Stokes operator initialized for selective use:
   div.initialize (system_mf_div_storage, std::vector (1,1), 
std::vector (1,0)); // div is StokesOperator: 
MatrixFreeOperators::Base 
for multiplication:
 BlockVectorType tb1(1),tb2(1);
 tb1.block(0)=system_rhs.block(0);
 tb2.block(0)=system_rhs.block(1);
 div.vmult(tb2, tb1);
Results in error  because vmult relays on local_apply that requires vector 
with 2 blocks.

2) I've created divergence operator (B from Stokes). However, vmult checks 
requires src and dst to have same size and throws an error. I'm not sure if 
pre and post processing of constraints will work on non-square matrix.

Michał Wichrowski

The Divergence operator is defined as follows:


template 
class DivergenceOperator : public MatrixFreeOperators::Base  >
{
public:
 typedef typename DoFHandler::active_cell_iterator CellIterator;
 typedef LinearAlgebra::distributed::Vector VectorType;

  DivergenceOperator():
 MatrixFreeOperators::Base()
 {};

 virtual void compute_diagonal()
   {}


private:
  void
  local_apply (const MatrixFree &data,
 VectorType  &dst,
   const VectorType&src,
   const std::pair &cell_range) 
const;


  virtual void apply_add(VectorType &dst,
 const VectorType &src) const;

};


template 
void
DivergenceOperator
::apply_add (VectorType   &dst,
 const VectorType &src) const
{
this->data->cell_loop (&DivergenceOperator::local_apply, this, dst, 
src);
}

template 
void
DivergenceOperator
::local_apply(const MatrixFree &data,
VectorType &dst,
const VectorType &src,
const std::pair &cell_range) const
 {

 typedef VectorizedArray vector_t;
 FEEvaluation velocity (data, 0);
 FEEvaluation pressure (data, 1);

 for (unsigned int cell=cell_range.first; cellhttp://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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Block problem +MatrixFree

[Michał Wichrowski]

W dniu niedziela, 1 października 2017 21:20:06 UTC+2 użytkownik Daniel 
Arndt napisał:
>
> Michael,
>
> I implemented something similar to what you proposed. You might want to 
> have a look at https://github.com/dealii/dealii/pull/5175
>
> Best,
> Daniel
>
> Thanks,  for now I'm using mine version (that seems to be working), I'll 
switch later to this one.

I've found another problem with off-diagonal blocks. The multiplication by 
B and BT is done by cell_loop  and it triggers assertion in when used 
inside multigrid. The problem is that my divergence operator did 
had adjust_ghost_range_if_necessary method. I resolved the problem by 
copying from Base operator, but I think that it may be added to MatrixFree 
class so that it had not have to be copied every time.  

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Matrix-Free framework for vector problems: no assertion for components

[Michał Wichrowski]

Dear deal.II developers,
I've found out that in Matrix-free framework there is no assertion that 
checks if number of components of finite element matches number components 
used in cell operations. By mistake I used scalar FE instead vector and it 
took me quite a long time to figure it out. The program worked and even 
covered in reasonable number of iterations. I think FEEvaluation is right 
place to put this assertion.

Michał Wichrowski

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Bug in FEEvaluationAccess::get_symmetric_gradient with Number =float

[Michał Wichrowski]

Dear all,
In following line:
line 4545 of /include/deal.II/matrix_free/fe_evaluation.h:
  VectorizedArray half = make_vectorized_array (0.5);

0.5 is interpreted as double and it leads to using 
make_vectorized_array. That is OK if Number=double but for any 
other type (eg. float) it results is compilation error. It should be 
replaced by:
  VectorizedArray half = make_vectorized_array (0.5);

Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Matrix-Free framework for vector problems: no assertion for components

[Michał Wichrowski]

W dniu niedziela, 8 października 2017 12:55:02 UTC+2 użytkownik Martin 
Kronbichler napisał:
>
> Dear Michal,
>
> Could you please be a bit more specific? What was the error that you 
> obtained? What vector did you pass into 
> FEEvaluation::read_dof_values/distribute_local_to_global? I guess that you 
> handed in a simple vector, e.g. LinearAlgebra::distributed::Vector and got 
> a segmentation fault. But to add the assertion I need more details. The 
> issues that in principle, we do support using multiple components out of a 
> scalar FE - but in that case the user must hand in a BlockVector or 
> std::vector. We could catch that case in the machinery, though.
>
> Thanks!
>
> Best,
> Martin
>
Martin,
There was no segmentation fault, program runned without any errors. I 
realised that something is wrong when I outputed results.  This is mine 
local apply,  VectorType =LinearAlgebra::distributed::Vector, problem 
(Stokes) consists of dim components of velocity and 1 component of 
pressure. 
Instead of using FESystem(FE_Q, dim) I used FE_Q for velocity.
After fixing this it work OK, so it is just a problem with assertion. I'll 
try to make minimal code reproducing error  (or rather lack of error).
  
template 
void
SymGradMass
::local_apply (const MatrixFree &data,
 VectorType  &dst,
 const VectorType&src,
 const std::pair &cell_range) const
{

  typedef VectorizedArray vector_t;
  FEEvaluation velocity (data, 0);  
//only velocity is used here 0-th component.

  for (unsigned int cell=cell_range.first; cellhttp://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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Matrix-free: Computing coefficient involving material ID

[Michał Wichrowski]

Dear all,
How can I compute coefficient based on material ID? I tried:

template
template
inline void SymGradMass
::evaluate_coefficient(
std::map& rho_map,
std::map& miu_map) {
return;
  const unsigned int n_cells = this->data->n_macro_cells();
  rho.resize  (n_cells);
  mu.resize  (n_cells);
  for (unsigned int cell=0; cell::n_array_elements ; ++element){
  const unsigned int  ID=this->data->get_cell_iterator(cell,element,0 
)->material_id();
  rho[cell][element] = rho_map.at(ID) ;
  mu[cell][element]  = miu_map.at(ID) ;
  }
}
}
But it triggers assertion in get_cell_iterator(..). If I understand it 
correctly, problem is that total number of cells is not divisible by 
VectorizedArray::n_array_elements and for one of  macro cells 
calling get_cell_iterator(cell,element,0 ) with element>1 makes no sense. 
However I have no idea how to avoid it.

Thanks for any help,
Michał

Error message:


An error occurred in line <1302> of file 
 in 
function
typename dealii::DoFHandler::cell_iterator dealii::MatrixFree::get_cell_iterator(unsigned int, unsigned int, unsigned int) const 
[with int dim = 2; Number = float; typename 
dealii::DoFHandler::cell_iterator = 
dealii::TriaIterator, 
false> >]
The violated condition was: 
(vector_number) < (irreg_filled)
Additional information: 
Index 1 is not in the half-open range [0,1).

Stacktrace:
---
#0  ./StokesMatrixFree: dealii::MatrixFree<2, 
float>::get_cell_iterator(unsigned int, unsigned int, unsigned int) const
#1  ./StokesMatrixFree: void SymGradMass<2, 2, 
float>::evaluate_coefficient(std::map, std::allocator > >&, std::map, 
std::allocator > >&)

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Bug in FEEvaluationAccess::get_symmetric_gradient with Number =float

[Michał Wichrowski]

I've tried to create pull request, but I've got error message:

mwichro@Preludio:~/lib/dealii$ git push origin  fix_feevaluation_float_inst
Username for 'https://github.com': mwichro
Password for 'https://mwic...@github.com': 
remote: Permission to dealii/dealii.git denied to mwichro.
fatal: unable to access 'https://github.com/dealii/dealii.git/': The 
requested URL returned error: 403


W dniu niedziela, 8 października 2017 12:58:50 UTC+2 użytkownik Martin 
Kronbichler napisał:
>
> Dear Michal,
>
> You are right, this is a bug and your fix is indeed correct. Would you 
> like to create a pull request fixing it? Here are some basic instructions: 
> http://dealii.org/participate.html including a video lecture. This way 
> you get the proper credit for it.
>
> Best,
> Martin
>
>

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Michał Wichrowski]

That's very bad news for me. This workaround will work for Laplace problem 
but not for my case. I'm dealing with Stokes problem and my method requires 
that second block of right hand side is zero, otherwise MinRes crashes and 
GMRES have to be used (that is far less effective). 

Can You give me some details what is the problem? 

W dniu piątek, 27 października 2017 15:02:29 UTC+2 użytkownik Wolfgang 
Bangerth napisał:
>
> On 10/27/2017 06:44 AM, Bruno Turcksin wrote: 
> > 
> > Could you add this explanation to the documentation and/or to a 
> tutorial. I 
> > think mentioning it step-37 would be great. 
>
> Yes, I was actually thinking the same. Maybe in "Possibilities for 
> extensions" 
> in the results section of step-37, like we have for other programs? 
>
> Cheers 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: 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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Michał Wichrowski]

I think that the problem can be avoided by splitting operator A into two 
parts: dirichlet nodes and other. The structure will be:
A1 0
0I
and the operator A1 applies matrix-vector product using fixed values on 
Dirichlet nodes.

If I got it right, the problem is with this part of code:
Line:  1186 of  operatrors.h:
// set zero Dirichlet values on the input vector (and remember the src 
and
// dst values because we need to reset them at the end)
for (unsigned int j=0; j
  
(subblock(src,j).local_element(edge_constrained_indices[j][i]),
  
 subblock(dst,j).local_element(edge_constrained_indices[j][i]));
subblock(const_cast(src),j).local_element(edge_constrained_indices[j][i]) = 0.;
  }
  }

By replacing 0 with Dirichlet boundary value and then applying the operator 
without Dirichlet constrains the multiplication by A1 can be done. 
Multiplication by identity for Dirichlet nodes is already implemented. 
This will require some changes in matrix-free framework. MatrixFree object 
will handle only homogeneous constrains and Operator will handle the 
non-homogeneous one.  The same change have to be done 
in vmult_interface_down  and  vmult_interface_up.


-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Michał Wichrowski]

W dniu sobota, 28 października 2017 18:52:57 UTC+2 użytkownik Martin 
Kronbichler napisał:
>
> Dear Michal, 
>
> You mean because once you apply an inhomogeneous Dirichlet condition on 
> the velocity you also get a contribution because the divergence is not 
> zero? 

 Exactly.

> You could work around that by applying a function in the whole 
> domain that satisfies the Dirichlet constraints and is divergence free, 
> in which case you should not get a contribution to right right hand side 
> of the divergence equation but only the momentum equation. On the other 
> hand, I'm not sure how easy it would be to construct such a function... 
>

Constructing such a function is not an option (for now), it is too 
expensive. 

One question about that way of imposing Dirichlet boundary conditions. I've 
modified step-37 following the same idea, but in different way:

  template 
  void LaplaceProblem::solve ()
  {
  SystemMatrixType system_matrix_no_dirichlet;
{
  typename MatrixFree::AdditionalData additional_data;
  additional_data.tasks_parallel_scheme =
MatrixFree::AdditionalData::none;
  additional_data.mapping_update_flags = (update_gradients | 
update_JxW_values |
  update_quadrature_points);
  std::shared_ptr >
  system_mf_storage(new MatrixFree());
  system_mf_storage->reinit (dof_handler, 
constraints_without_dirichlet, QGauss<1>(fe.degree+1),
 additional_data);
  system_matrix_no_dirichlet.initialize (system_mf_storage);
}
system_matrix_no_dirichlet.evaluate_coefficient(Coefficient());

//
constraints.distribute(solution);
LinearAlgebra::distributed::Vector residual=system_rhs;
system_matrix_no_dirichlet.vmult(residual, solution);
system_rhs -= residual;
LinearAlgebra::distributed::Vector solution_update=system_rhs;
solution_update=0;
..

cg.solve (system_matrix, solution_update, system_rhs,
  preconditioner);
solution+= solution_update;

constraints.distribute(solution);


In this way I don't have to modify my operator, bum I'm not sure if it will 
work for all cases. For now it works and the results looks OK.
 

>
> Regarding the workaround: It could work by the modification you 
> suggested, but we cannot really fix that in the library because the 
> operators are meant for linear operators and not affine ones. The 
> problematic issue is that there are some cases where you want to apply 
> an operator with a non-zero constraint and others where you don't. For 
> example, all Krylov solvers expect you to provide a matrix-vector for 
> the linear operator, not the affine one as you would get when you have 
> inhomogeneous constraints inside the mat-vec. In particular, I do not 
> really understand the part in your second post where you talk about the 
> vmult_interface_{up,down} functions because those are only used for 
> multigrid where the library always assumes to work on homogeneous 
> problems. 
>
> You're right,  forget about what I written previously. I've totally 
misunderstood a few things.
 

> I have not gone down this route and as far as I know, nobody else has 
> previously - including the matrix-based solvers we have available in 
> deal.II whose way of operation resembles the setup I described in my 
> previous post, but do it inside the 
> ConstraintMatrix::distribute_local_to_global call - so I cannot say how 
> much work that would be. We would be happy to provide a solution that is 
> generic if you find one, but I'm not sure I will be able to help much 
> along this direction. 
>
> Best, 
> Martin 
>

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Non-homogeneus boundary conditions with-matrix-free not working

[Michał Wichrowski]

Dear Martin,
I finally managed to make it work by following Your idea. Thanks!

Michał 

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] LinearOperator -- inverse_operator : Payload for LinearAlgebra::distributed::Vector

[Michał Wichrowski]

Dear all, 
I'm trying to use iterative inverse with distributed deal.II vector. From 
following lines 

const auto S= linear_operator (mg_sc_matrices[level]);
const auto Shat = 
inverse_operator(S,solver,SCPreconditionMGs[level]);

 I get compilation errors:


/home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
note: candidate: template dealii::LinearOperator dealii::inverse_operator(const 
dealii::LinearOperator&, Solver&, const 
Preconditioner&)
 inverse_operator(const LinearOperator &op,
 ^
/home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
note:   template argument deduction/substitution failed:
/home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
note:   mismatched types 
‘dealii::LinearAlgebra::distributed::Vector’ and 
‘dealii::internal::LinearOperator::EmptyPayload’
const auto Shat = 
inverse_operator(S,solver,SCPreconditionMGs[level]);
   ^
/home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
note:   ‘const 
dealii::LinearOperator, 
dealii::LinearAlgebra::distributed::Vector, 
dealii::internal::LinearOperator::EmptyPayload>’ is not derived from ‘const 
dealii::LinearOperator >’
CMakeFiles/StokesMatrixFree.dir/build.make:62: recipe for target 
'CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o' failed
make[2]: *** [CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o] Error 1


I think that proper Payload is needed, but I didn't found right one.

Michał 

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] LinearOperator -- inverse_operator : Payload for LinearAlgebra::distributed::Vector

[Michał Wichrowski]

Sorry, I totally forgotten about types:

typedef LinearAlgebra::distributed::Vector LevelVectorType;
SolverCG solver(solver_control);
std::vector  >  > SCPreconditionMGs;



W dniu sobota, 17 lutego 2018 06:32:07 UTC+1 użytkownik Matthias Maier 
napisał:
>
> Hi Michał, 
>
> It would be helpful if you could sent us a small, minimal example that 
> shows the problem (in particular, what type is "LevelVectorType" and 
> "SCPreconditoinMGs", "solver"). 
> The example should be as short as possible and doesn't need to be 
> functional / or do anything at all - we essentially only need to be able 
> to figure out the types you try to wrap into the linear_operator object. 
>
> Best, 
> Matthias 
>
>
> On Tue, Feb 13, 2018, at 12:47 CST, Michał Wichrowski  > wrote: 
>
> > Dear all, 
> > I'm trying to use iterative inverse with distributed deal.II vector. 
> From 
> > following lines 
> > 
> > const auto S= linear_operator 
> (mg_sc_matrices[level]); 
> > const auto Shat = 
> > inverse_operator(S,solver,SCPreconditionMGs[level]); 
> > 
> >  I get compilation errors: 
> > 
> > 
> > /home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
> > note: candidate: template > Preconditioner, class Range, class Domain> 
> dealii::LinearOperator > Range, Payload> dealii::inverse_operator(const 
> > dealii::LinearOperator&, Solver&, const 
> > Preconditioner&) 
> >  inverse_operator(const LinearOperator &op, 
> >  ^ 
> > /home/mwichro/lib/deal.II/include/deal.II/lac/linear_operator.h:678:1: 
> > note:   template argument deduction/substitution failed: 
> > 
> /home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
> > note:   mismatched types 
> > ‘dealii::LinearAlgebra::distributed::Vector’ and 
> > ‘dealii::internal::LinearOperator::EmptyPayload’ 
> > const auto Shat = 
> > inverse_operator(S,solver,SCPreconditionMGs[level]); 
> >^ 
> > 
> /home/mwichro/deal-ii-projects/StokeMatrixFree/StokesMatrixFree.cc:1051:55: 
> > note:   ‘const 
> > 
> dealii::LinearOperator, 
> > dealii::LinearAlgebra::distributed::Vector, 
> > dealii::internal::LinearOperator::EmptyPayload>’ is not derived from 
> ‘const 
> > dealii::LinearOperator > dealii::LinearAlgebra::distributed::Vector >’ 
> > CMakeFiles/StokesMatrixFree.dir/build.make:62: recipe for target 
> > 'CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o' failed 
> > make[2]: *** [CMakeFiles/StokesMatrixFree.dir/StokesMatrixFree.cc.o] 
> Error 1 
> > 
> > 
> > I think that proper Payload is needed, but I didn't found right one. 
> > 
> > Michał 
>

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] source/base/mpi.cc headers

[Michał Wichrowski]

Dear all,
I had some trouble while compiling deal.II on my cluster:

[ 70%] Building CXX object 
source/fe/CMakeFiles/obj_fe_release.dir/fe_tools.cc.o
/home/mwichro/lib/dealii-9.0.0/source/base/mpi.cc: In function 'int 
dealii::Utilities::MPI::create_group(ompi_communicator_t* const&, 
ompi_group_t* const&, int, ompi_communicator_t**)':
/home/mwichro/lib/dealii-9.0.0/source/base/mpi.cc:129:7: error: 'iota' is 
not a member of 'std'
   std::iota(grp_pids.begin(), grp_pids.end(), 0);
   ^~~
/home/mwichro/lib/dealii-9.0.0/source/base/mpi.cc: In function 'int 
dealii::Utilities::MPI::create_group(ompi_communicator_t* const&, 
ompi_group_t* const&, int, ompi_communicator_t**)':
/home/mwichro/lib/dealii-9.0.0/source/base/mpi.cc:129:7: error: 'iota' is 
not a member of 'std'
   std::iota(grp_pids.begin(), grp_pids.end(), 0);
   ^~~

I fixed that by adding 
 #include  
to  source/base/mpi.cc, according 
to https://cdcvs.fnal.gov/redmine/projects/larsoft/wiki/Updating_code_to_GCC_6

The problem will probably also  appear in developement version of deal.II. 


Best, 
Michal


-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Compilation problem with Development sources

[Michał Wichrowski]

Dear all,
I had a problem while compiling deal.II source code from the git repository:

[ 72%] Building CXX object 
source/algorithms/CMakeFiles/obj_algorithms_debug.dir/operator.cc.o
In file included from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/any_data.h:21:0,
 from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/newton.h:20,
 from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/newton.templates.h:20,
 from 
/home/mwichro/lib/dealii/source/algorithms/operator.cc:17:
/home/mwichro/lib/dealii/include/deal.II/base/utilities.h: In function 
‘constexpr unsigned int dealii::Utilities::pow(unsigned int, int)’:
/home/mwichro/lib/dealii/include/deal.II/base/exceptions.h:1242:72: error: 
call to non-constexpr function ‘void 
dealii::deal_II_exceptions::internals::issue_error_noreturn(dealii::deal_II_exceptions::internals::ExceptionHandling,
 
const char*, int, const char*, const char*, const char*, ExceptionType) 
[with ExceptionType = dealii::StandardExceptions::ExcMessage]’
   ::dealii::deal_II_exceptions::internals::issue_error_noreturn( \
^
/home/mwichro/lib/dealii/include/deal.II/base/utilities.h:357:5: note: in 
expansion of macro ‘Assert’
 Assert(iexp >= 0, ExcMessage("The exponent must not be negative!"));
 ^
source/algorithms/CMakeFiles/obj_algorithms_debug.dir/build.make:54: recipe 
for target 
'source/algorithms/CMakeFiles/obj_algorithms_debug.dir/operator.cc.o' failed
make[2]: *** 
[source/algorithms/CMakeFiles/obj_algorithms_debug.dir/operator.cc.o] Error 
1
CMakeFiles/Makefile2:4090: recipe for target 
'source/algorithms/CMakeFiles/obj_algorithms_debug.dir/all' failed
make[1]: *** [source/algorithms/CMakeFiles/obj_algorithms_debug.dir/all] 
Error 2
Makefile:116: recipe for target 'all' failed
make: *** [all] Error 2



Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Compilation problem with Development sources

[Michał Wichrowski]

I've just double checked,  I'm using the latest version:

mwichro@major:~/lib/dealii$ git checkout 
Your branch is up-to-date with 'origin/master'.
mwichro@major:~/lib/dealii$ make
[  0%] Built target expand_instantiations_exe
[  0%] Built target doxygen_headers
[  3%] Built target obj_boost_serialization_debug
[  5%] Built target obj_boost_serialization_release
[  5%] Built target obj_boost_system_debug
[  5%] Built target obj_boost_system_release
[  7%] Built target obj_tbb_debug
[  9%] Built target obj_tbb_release
[  9%] Built target obj_umfpack_DI_ASSEMBLE_debug
[ 10%] Built target obj_umfpack_DI_ASSEMBLE_release
[ 10%] Built target obj_umfpack_DI_SOLVE_debug
[ 10%] Built target obj_umfpack_DI_SOLVE_release
[ 10%] Built target obj_umfpack_DI_STORE_debug
[ 10%] Built target obj_umfpack_DI_STORE_release
[ 10%] Built target obj_umfpack_DI_TRIPLET_MAP_NOX_debug
[ 10%] Built target obj_umfpack_DI_TRIPLET_MAP_NOX_release
[ 10%] Built target obj_umfpack_DI_TRIPLET_MAP_X_debug
[ 10%] Built target obj_umfpack_DI_TRIPLET_MAP_X_release
[ 10%] Built target obj_umfpack_DI_TRIPLET_NOMAP_X_debug
[ 10%] Built target obj_umfpack_DI_TRIPLET_NOMAP_X_release
[ 10%] Built target obj_umfpack_DI_TSOLVE_debug
[ 10%] Built target obj_umfpack_DI_TSOLVE_release
[ 11%] Built target obj_umfpack_DL_ASSEMBLE_debug
[ 11%] Built target obj_umfpack_DL_ASSEMBLE_release
[ 11%] Built target obj_umfpack_DL_SOLVE_debug
[ 11%] Built target obj_umfpack_DL_SOLVE_release
[ 11%] Built target obj_umfpack_DL_STORE_debug
[ 11%] Built target obj_umfpack_DL_STORE_release
[ 11%] Built target obj_umfpack_DL_TRIPLET_MAP_NOX_debug
[ 11%] Built target obj_umfpack_DL_TRIPLET_MAP_NOX_release
[ 11%] Built target obj_umfpack_DL_TRIPLET_MAP_X_debug
[ 11%] Built target obj_umfpack_DL_TRIPLET_MAP_X_release
[ 11%] Built target obj_umfpack_DL_TRIPLET_NOMAP_X_debug
[ 11%] Built target obj_umfpack_DL_TRIPLET_NOMAP_X_release
[ 11%] Built target obj_umfpack_DL_TSOLVE_debug
[ 13%] Built target obj_umfpack_DL_TSOLVE_release
[ 13%] Built target obj_umfpack_GENERIC_debug
[ 13%] Built target obj_umfpack_GENERIC_release
[ 15%] Built target obj_umfpack_I_UMFPACK_debug
[ 17%] Built target obj_umfpack_I_UMFPACK_release
[ 20%] Built target obj_umfpack_I_UMF_debug
[ 23%] Built target obj_umfpack_I_UMF_release
[ 26%] Built target obj_umfpack_L_UMFPACK_debug
[ 28%] Built target obj_umfpack_L_UMFPACK_release
[ 31%] Built target obj_umfpack_L_UMF_debug
[ 34%] Built target obj_umfpack_L_UMF_release
[ 35%] Built target obj_amd_global_debug
[ 35%] Built target obj_amd_global_release
[ 35%] Built target obj_amd_int_debug
[ 36%] Built target obj_amd_int_release
[ 38%] Built target obj_amd_long_debug
[ 39%] Built target obj_amd_long_release
[ 40%] Built target obj_muparser_debug
[ 41%] Built target obj_muparser_release
[ 41%] Built target obj_matrix_free_inst
[ 41%] Built target obj_matrix_free_release
[ 41%] Built target obj_meshworker_inst
[ 42%] Built target obj_meshworker_release
[ 42%] Built target obj_opencascade_inst
[ 42%] Built target obj_opencascade_release
[ 45%] Built target obj_fe_inst
[ 48%] Built target obj_fe_release
[ 48%] Built target obj_distributed_inst
[ 50%] Built target obj_distributed_release
[ 50%] Built target obj_physics_elasticity_inst
[ 50%] Built target obj_physics_elasticity_release
[ 50%] Built target obj_sundials_inst
[ 50%] Built target obj_sundials_release
[ 50%] Built target obj_physics_inst
[ 51%] Built target obj_physics_release
[ 51%] Built target obj_non_matching_inst
[ 51%] Built target obj_non_matching_release
[ 52%] Built target obj_base_inst
[ 56%] Built target obj_base_release
[ 56%] Built target obj_algorithms_inst
[ 56%] Built target obj_algorithms_release
[ 58%] Built target obj_lac_inst
[ 60%] Built target obj_lac_release
[ 61%] Built target obj_dofs_inst
[ 63%] Built target obj_dofs_release
[ 65%] Built target obj_numerics_inst
[ 67%] Built target obj_numerics_release
[ 68%] Built target obj_multigrid_inst
[ 68%] Built target obj_multigrid_release
[ 69%] Built target obj_particle_inst
[ 69%] Built target obj_particle_release
[ 69%] Built target obj_hp_inst
[ 70%] Built target obj_hp_release
[ 71%] Built target obj_grid_inst
[ 72%] Built target obj_grid_release
[ 72%] Built target obj_gmsh_inst
[ 72%] Built target obj_gmsh_release
[ 72%] Built target obj_differentiation_ad_inst
[ 72%] Built target obj_differentiation_ad_release
[ 72%] Built target deal_II
[ 72%] Building CXX object 
source/algorithms/CMakeFiles/obj_algorithms_debug.dir/operator.cc.o
In file included from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/any_data.h:21:0,
 from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/newton.h:20,
 from 
/home/mwichro/lib/dealii/include/deal.II/algorithms/newton.templates.h:20,
 from 
/home/mwichro/lib/dealii/source/algorithms/operator.cc:17:
/home/mwichro/lib/dealii/include/deal.II/base/utilities.h: In function 
‘constexpr unsigned int dealii::Utilities::pow(unsigned int, int)’:

Re: [deal.II] Compilation problem with Development sources

[Michał Wichrowski]

W dniu poniedziałek, 17 września 2018 16:10:32 UTC+2 użytkownik Jean-Paul 
Pelteret napisał:
>
> Dear Michał,
>
> I think that calling *git checkout* is insufficient. You probably want to 
> call* git pull *(or* git fetch && git rebase origin/master*) to fetch and 
> merge the most current state of the remote repository.
>
> Best,
> Jean-Paul
>
> Thanks, the compilation seems to get though that problems. I'm still 
waiting for compiler to finish but I think it will work.
Michał 

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Compilation error: use of deleted function

[Michał Wichrowski]

Dear all,

I have errors when compiling my code with the devellpment version of 
deal.II .  The same problems occurs also with 9.0.0, the code works well 
with 9.0.0-pre. I  have no idea how to resolve these problems.

Michał 

/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:137:7: 
error: use of deleted function ‘dealii::mg::Signals::Signals(const 
dealii::mg::Signals&)’
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
note: ‘dealii::mg::Signals::Signals(const dealii::mg::Signals&)’ is 
implicitly deleted because the default definition would be ill-formed:
   struct Signals
  ^
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function 
‘boost::signals2::signal::signal(const 
boost::signals2::signal&)’
In file included from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/slot.hpp:21:0,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/connection.hpp:24,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/signal.hpp:22,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2.hpp:19,
 from 
/home/mwichro/lib/deal.II/include/deal.II/lac/solver.h:27,
 from 
/home/mwichro/lib/deal.II/include/deal.II/lac/solver_cg.h:26,
 from 
/home/mwichro/deal.II-projects/StokeMatrixFree/StokesMatrixFree.cc:8:
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/detail/signal_template.hpp:599:11:
 
note: ‘boost::signals2::signal::signal(const 
boost::signals2::signal&)’ is implicitly declared as deleted 
because ‘boost::signals2::signal’ declares a move constructor 
or move assignment operator
 class BOOST_SIGNALS2_SIGNAL_CLASS_NAME(BOOST_SIGNALS2_NUM_ARGS)
   ^
In file included from 
/home/mwichro/deal.II-projects/StokeMatrixFree/StokesMatrixFree.cc:34:0:
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function 
‘boost::signals2::signal::signal(const 
boost::signals2::signal&)’
   struct Signals
  ^
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function ‘boost::signals2::signal::signal(const boost::signals2::signal&)’
In file included from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/slot.hpp:21:0,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/connection.hpp:24,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/signal.hpp:22,
 from 
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2.hpp:19,
 from 
/home/mwichro/lib/deal.II/include/deal.II/lac/solver.h:27,
 from 
/home/mwichro/lib/deal.II/include/deal.II/lac/solver_cg.h:26,
 from 
/home/mwichro/deal.II-projects/StokeMatrixFree/StokesMatrixFree.cc:8:
/home/mwichro/lib/deal.II/include/deal.II/bundled/boost/signals2/detail/signal_template.hpp:599:11:
 
note: ‘boost::signals2::signal::signal(const 
boost::signals2::signal&)’ is implicitly declared 
as deleted because ‘boost::signals2::signal’ 
declares a move constructor or move assignment operator
 class BOOST_SIGNALS2_SIGNAL_CLASS_NAME(BOOST_SIGNALS2_NUM_ARGS)
   ^
In file included from 
/home/mwichro/deal.II-projects/StokeMatrixFree/StokesMatrixFree.cc:34:0:
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function ‘boost::signals2::signal::signal(const boost::signals2::signal&)’
   struct Signals
  ^
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function ‘boost::signals2::signal::signal(const boost::signals2::signal&)’
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function ‘boost::signals2::signal::signal(const boost::signals2::signal&)’
/home/mwichro/lib/deal.II/include/deal.II/multigrid/multigrid.h:53:10: 
error: use of deleted function ‘boost::signals2::signal::signal(const boost::signals2::signal&)’

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compilation error: use of deleted function

[Michał Wichrowski]

>
> Michal, 
> the error is in your code, but since I don't have it, I can't tell 
> where. I suspect that you are trying to copy an object of type 
> Multigrid, in StokeMatrixFree/StokesMatrixFree.cc:34 of your project. 
> This copy may happen implicitly if you call a function that takes a 
> Multigrid object by value, rather than by reference. (Or does so for an 
> object that has a Multigrid object as a member variable.) 
>
> Best 
>   W. 
>

Yep, I am trying to copy Multigrid object on purpose, I found the line 
corresponding to the  problem:
Line 762:
  std::vector >  
multigrids_sc(triangulation.n_global_levels(),
  mg_sc);

 I need a vector of objects of type Multigrid, how can I fix this?

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compilation error: use of deleted function

[Michał Wichrowski]

W dniu wtorek, 25 września 2018 15:30:06 UTC+2 użytkownik Wolfgang Bangerth 
napisał:
>
> On 09/25/2018 05:56 AM, Michał Wichrowski wrote: 
> > 
> > Yep, I am trying to copy Multigrid object on purpose, I found the line 
> > corresponding to the  problem: 
> > Line 762: 
> >std::vector >   
> > multigrids_sc(triangulation.n_global_levels(), 
> >mg_sc); 
> > 
> >   I need a vector of objects of type Multigrid, how can I fix this? 
>
> Do you actually need to copy the Multigrid objects, or do you just want to 
> store them in a vector? (std::vector of course requires its objects to be 
> copiable.) 
>
> If you don't actually need copies, then just use a 
> std::vector>. 
>
> Best 
>   W. 
>

Yes, I need triangulation.n_global_levels() copies of multigrid. I modify 
each of them in next line:

  std::vector >  
multigrids_sc(triangulation.n_global_levels(),
  mg_sc);
  for (unsigned int level = 0; levelhttp://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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compilation error: use of deleted function

[Michał Wichrowski]

>
> On 09/25/2018 07:38 AM, Michał Wichrowski wrote: 
> > 
> > 
> > Yes, I need triangulation.n_global_levels() copies of multigrid. I 
> modify each 
> > of them in next line: 
> > 
> >std::vector >   
> > multigrids_sc(triangulation.n_global_levels(), 
> >mg_sc); 
> >for (unsigned int level = 0; level ++level) 
> >multigrids_sc[level].reinit(0, level); 
>
> I know you need this many objects, but what I meant is: do you actually 
> need 
> to copy one object to another? If you use a 
> std::vector>, 
> you can also have as many copies as you showed. 
>
> Best 
>   W. 
>
> Ok, now I understand the idea. I do not need to copy a specific object, 
just several Multigrid objects, so it should work.
If I understand correctly, I should do something like this:
  std::vector  > >
multigrids_sc_pointers(triangulation.n_global_levels());
  for (unsigned int level = 0; level(
mg_matrix_sc,
mg_coarse_sc,
mg_transfer_sc,
mg_smoother_sc,
mg_smoother_sc);
It does not work, I have almost no experience with unique pointers and I 
have no idea how to do it. Could you help me? 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Compilation error: use of deleted function

[Michał Wichrowski]

W dniu wtorek, 25 września 2018 20:59:41 UTC+2 użytkownik Wolfgang Bangerth 
napisał:
>
> On 09/25/2018 08:33 AM, Michał Wichrowski wrote: 
> >multigrids_sc_pointers[level] =  new Multigrid( 
> >  mg_matrix_sc, 
> >  mg_coarse_sc, 
> > mg_transfer_sc, 
> > mg_smoother_sc, 
> > mg_smoother_sc); 
>
> This won't compile. You'll have to write this as 
>
> multigrids_sc_pointers[level] = 
>std_cxx14::make_unique>( 
> mg_matrix_sc, 
> mg_coarse_sc, 
> mg_transfer_sc, 
> mg_smoother_sc, 
> mg_smoother_sc); 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

It works
Thank You
Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Computing min and max eigenvalues

[Michał Wichrowski]

Dear all,
I need to compute the minimum and maximum eigenvalues of the matrix-free 
operator (vmult is only available ). For the largest eigenvalue, I'm using 
the power method, but I have some problems with minimal eigenvalue. The 
condition number of the operator is close to 1, so I tried using the power 
method on inverse obtained by CG. This solution is time-consuming and for 
some cases I got strange results.
I have seen that CG can output eigenvalues: 
  solver_outer.connect_eigenvalues_slot( [](const std::vector<  double > & 
vec)->void { std::cout<< "  "

[deal.II] Re: Computing min and max eigenvalues

[Michał Wichrowski]

W dniu czwartek, 18 października 2018 17:59:40 UTC+2 użytkownik Bruno 
Turcksin napisał:
>
> Michal,
>
> Have you tried to use ARPACK to get the smallest eigenvalues? That's what 
> we are doing in one of our code and it works pretty well. We have a patch 
> here 
> 
>  
> so that you don't need to provide the inverse of the operator to use 
> ARPACK. You still need to provide a mass matrix but it's no used.
>
> Best,
>
> Bruno
>

Thanks a lot!
In the meantime, I found that  I was probably using the wrong solver to get 
inverse (the system may be not SPD) using GMRes resulted in a 
significant speedup. However, I still have strange results, including the 
smallest eigenvalue being larger than the largest one. I think using ARPACK 
instead of my own code will solve the problem, thank you.
Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Computing min and max eigenvalues

[Michał Wichrowski]

Dear Bruno,
I've used ARPACK and the computation time has dropped signifficantly. But I 
have some strange results, I have following code:

  IterationNumberControl arpack_control(30, 1e-8);

  typename ArpackSolver::AdditionalData arpack_data_smallest(30,
   ArpackSolver:: 
smallest_magnitude,
   false);
  typename ArpackSolver::AdditionalData arpack_data_largest(30,
   ArpackSolver:: 
largest_magnitude,
   false);
ArpackSolver arpack_smallest (arpack_control,arpack_data_smallest);
arpack_smallest.solve(shat,mass, s_hatS_inverse, eigenvalues, 
eigenvectors_small, 1 );

pcout<<"\t Aparck small "<<"coverged in "
< eigenvectors_large(2, rhs);
ArpackSolver arpack_largest (arpack_control,arpack_data_largest);
arpack_largest.solve(shat,mass, s_hatS_inverse, eigenvalues, 
eigenvectors_large, 1 );

pcout<<"\t Aparck large "<<"coverged in "
<::AdditionalData data(0.);
  EigenPower eigen(solver_control_eigen,mem,data );
  double l_max=1;
  double l_min=0;
//   rhs=1;  //starting from 1 gives same results
  rhs= eigenvectors_large[0];
  internal::PreconditionChebyshevImplementation::set_initial_guess(rhs);
  eigen.solve(l_max,shat, rhs );
  sol=1;

pcout<<"\t Eigen power  "<<"coverged in "

[deal.II] Re: Computing min and max eigenvalues

[Michał Wichrowski]

PS. 
I noticed that in some cases the "smallest" eigenvalue does not 
match maximum eigenvalue obtained from power method.

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Computing min and max eigenvalues

[Michał Wichrowski]

Dear Bruno

>
> > 2) The largest eigenvalue is computed instead of the largest one, and, I 
> guess, the smallest eigenvalue is computed instead of the largest one. 
> If you use ARPACK through deal without patching it, then the only 
> method available is shift-inverse. The thing is that when you know the 
> largest eigenvalue of the inverse you can get the smallest eigenvalue 
> of the operator. So yes, the names are reversed. 
>
Yup,  using the patch solved this issue.  At least the smallest eigenvalue 
is smaller than the largest.

>
> > 3) I get complex eigenvalues in some cases and I am quite sure, that all 
> eigenvalues of the system are real. 
> I've not had this problem. Maybe it's because of numerical error? 

Eigenvalue:  (1.03342,0.0069945) with arpack tolerance 1e-8. Is it possible?

>
> > I noticed that in some cases the "smallest" eigenvalue does not match 
> maximum eigenvalue obtained from power method. 
> Are you sure that the power method has converged? Also, try using 
> ARPACK's regular mode instead of shift-inverse. I think it might be 
> better to find the largest eigenvalue. 
>
 The deal.II power method converged with residual: 9.5052e-05 I switched to 
regular mode and obtained 1.04192 from arpack and 1.03981 from power 
method. The tolerance of arpack was set to 1e-8. The difference is growing 
with mesh size,  so I assume the tolerance is set for eigenvector, 
resulting in much smaller tolerance for the eigenvalue.

Also, I got the smallest eigenvalues equal to  0 in some cases, while I'm 
sure that the matrix is not singular. I can solve a linear system with the 
same matrix using  Krylov solvers without a problem.

Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Computing min and max eigenvalues

[Michał Wichrowski]

Earl,

I am posing this question purely out of personal interest: what do you mean 
> by maximum and minimum eigenvalue? I am confused by this.
>
 
>
The eigenvalues with largest and smallest magnitude. My problem is 
symmetric positive-definite, thus the spectrum is real and positive, so in 
my case "minimum eigenvalue" makes sense. 

Bruno,
I think I finally have a solution. By increasing Arnoldi space size, I 
obtained results that seems ok.

I think the documentation of arpack solver is not valid: 
https://www.dealii.org/developer/doxygen/deal.II/classArpackSolver.html#afdc3aa9d761c43b5b4132e731c6191a5
A The operator for which we want to compute eigenvalues. Actually, this 
parameter is entirely unused.*Actually, this parameter is entirely unused*.The 
matrix A is obviously used. 

Also, calling arpack with WhichEigenvalues::both_ends always results in a 
crash.
Thanks for help,
Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] copy_triangulation removes limit_level_difference_at_vertices flag

[Michał Wichrowski]

Dear all,
I noticed that if I do the following:

parallel::distributed::Triangulation  
triangulation(MPI_COMM_WORLD,
parallel::distributed::Triangulation::limit_level_difference_at_vertices,
  
 parallel::distributed::Triangulation::construct_multigrid_hierarchy);

  Triangulation  other_tria;
  //do something with other tria
triangulation.copy_triangulation(other_tria);

dof_handler_u.distribute_dofs(fe_velocity);
dof_handler_u.distribute_mg_dofs();
I got an error: 

An error occurred in line <1275> of file 
 in function
void dealii::DoFHandler::distribute_mg_dofs() [with int 
dim = 2; int spacedim = 2]
The violated condition was: 
((tria->get_mesh_smoothing() & Triangulation::limit_level_difference_at_vertices) != Triangulation::none)
Additional information: 
The mesh smoothing requirement 'limit_level_difference_at_vertices' has 
to be set for using multigrid!

This is fixed by   Triangulation  
other_tria(Triangulation::limit_level_difference_at_vertices);

Is this intendent to work like that or it is a bug?

Michał

-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: copy_triangulation removes limit_level_difference_at_vertices flag

[Michał Wichrowski]

Done: https://github.com/dealii/dealii/issues/7581

W dniu wtorek, 8 stycznia 2019 23:42:06 UTC+1 użytkownik Daniel Arndt 
napisał:
>
> Michal,
>
> not setting `limit_level_difference_at_vertices` when copying a 
> parallel::distributed::Triangulation can be considered an oversight/bug.
> Do you want to create a PR for fixing this?
>
> 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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Multigrid starting from other level than finest does not work with adaptive refinement

[Michał Wichrowski]

Dear all,
I am able to reproduce the problem with modified step-37, I expect the 
similar problem in case of multigrid with sparse matrices. The problem 
occurs only if mesh is refined adaptively.
I added  refine_mesh() function and  modified solve():

//
//work on level 3, that is not the finest level
unsigned int lvl=3;
  mg.reinit(0, lvl);
//..
//..
//Set new rhs vector of propper size
LinearAlgebra::distributed::Vector rhs_lvl;
mg_matrices[lvl].get_matrix_free()-> initialize_dof_vector(rhs_lvl);
LinearAlgebra::distributed::Vector sol=rhs_lvl;
//fill rhs_lvl with anything:
rhs_lvl=1;
  pcout  << " rhs size (1) :  "< of file 
 
in function
void 
dealii::MGLevelGlobalTransfer >::copy_to_mg(const dealii::DoFHandler&, 
dealii::MGLevelObject >&, const 
dealii::LinearAlgebra::distributed::Vector&, bool) const [with int 
dim = 3; Number2 = float; int spacedim = 3; Number = float]
The violated condition was: 
(ghosted_global_vector.local_size()) == (src.local_size())
Additional information: 
Dimension 27150 not equal to 4913.

Stacktrace:
---
#0  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.1.0-pre: void 
dealii::MGLevelGlobalTransfer >::copy_to_mg<3, float, 3>(dealii::DoFHandler<3, 
3> const&, 
dealii::MGLevelObject >&, 
dealii::LinearAlgebra::distributed::Vector const&, bool) const
#1  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.1.0-pre: void 
dealii::MGLevelGlobalTransfer >::copy_to_mg<3, float, 3>(dealii::DoFHandler<3, 
3> const&, 
dealii::MGLevelObject >&, 
dealii::LinearAlgebra::distributed::Vector const&) const
#2  ./step-37: void 
dealii::internal::PreconditionMGImplementation::vmult<3, 
dealii::LinearAlgebra::distributed::Vector, dealii::MGTransferMatrixFree<3, float>, 
dealii::LinearAlgebra::distributed::Vector >(std::vector const*, 
std::allocator const*> > const&, 
dealii::Multigrid >&, dealii::MGTransferMatrixFree<3, float> 
const&, dealii::LinearAlgebra::distributed::Vector&, 
dealii::LinearAlgebra::distributed::Vector const&, bool, dealii::mg::Signals const&, ...)
#3  ./step-37: void dealii::PreconditionMG<3, 
dealii::LinearAlgebra::distributed::Vector, dealii::MGTransferMatrixFree<3, float> 
>::vmult 
>(dealii::LinearAlgebra::distributed::Vector&, 
dealii::LinearAlgebra::distributed::Vector const&) const
#4  ./step-37: Step37::LaplaceProblem<3>::solve()
#5  ./step-37: Step37::LaplaceProblem<3>::run()
#6  ./step-37: main


-- 
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2009 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
 * 2009-2012, updated to MPI version with parallel vectors in 2016
 */


#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 
#include 


#include 
#include 


namespace Step37
{
  using namespace dealii;


  const unsigned int degree_finite_element = 2;
  const unsigned int dimension = 3;



  template 
  class Coefficient : public Function
  {
  public:
Coefficient()
  : Function()
{}

virtual double value(const Point & p,
 const unsigned int component = 0) const override;

template 
number value(const Point &p,
 const unsigned intcomponent = 0) const;
  };



  template 
  template 
  number Coefficient::value(const Point &p,
 const unsigned int /*component*/) const
  {
return 1. / (0.05 + 2. * p.square());
  }



  template 
  double Coefficient::value(const Point & p,
 const unsigned int component) const
  {
return value(p, component);
  }



  template 
  class LaplaceOperator
: public MatrixFreeOperators::
Base>
  {
  public:
using value_type = number;

LaplaceOperator();

 

[deal.II] DG Matrix-free: getting face iterator form

[Michał Wichrowski]

Dear all,
I'm trying to use matrix-free with discontinuous Galerkin elements for 
variable coefficient problem.  The face integral require access to value of 
coefficient. For the cell I've done the following thing:
//Vectors of values, stored in object:
  AlignedVector< VectorizedArray > mu;
  AlignedVector< VectorizedArray > rho;
//Filling the vectors:
  const unsigned int n_cells = this->data->n_macro_cells();
  rho.resize  (n_cells);
  mu.resize  (n_cells);
  for (unsigned int cell=0; celldata->n_components_filled(cell) 
; ++element){
  const unsigned int  ID=this->data->get_cell_iterator(cell,element,0 
)->material_id();
  rho[cell][element] = rho_map.at(ID) ;
  mu[cell][element]  = miu_map.at(ID) ;
  }
}
For the faces I am trying to do something like that:
//Vectors of values, stored in object:
  AlignedVector< VectorizedArray > mu_inner_faces;

//Filling the vectors:
  const unsigned int n_faces=this->data->n_inner_face_batches();
  mu_inner_faces.resize(n_faces);
  for (unsigned int face=0; face::n_array_elements > 
face2cell_info=this->data->get_face_info(face);

  for(unsigned int element=0;
  elementdata->n_active_entries_per_face_batch(face) ;
  ++element){
  const unsigned int interior_cell = face2cell_info.cells_interior[element];
const unsigned int exterior_cell = 
face2cell_info.cells_exterior[element];
   
   ??

  }
  }

If I understood correctly interior_cell and exterior_cell will be numbers 
of cell.
 I have no idea how to get cell iterators from that information. I even 
found that MAtrixFree class contains attribute
std::vector< std::pair< unsigned int, unsigned int > >  cell_level_index 

but unfortunately it is private and  there is no  getter for that. 

I will also need to get user_index from each face and create from that 
another AlignedVector of parameters. This operations will be done only once 
so I do not care much about efficiency.
Best, 
Michal

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b3c0fccd-7a76-4d90-a573-74719d281c37%40googlegroups.com.

Re: [deal.II] DG Matrix-free: getting face iterator form

[Michał Wichrowski]

> Dear Martin
>
> We do not offer convenience access for all data structures on faces yet, 
> so feel free to suggest what you need.
>
1)   functions similar to get_cell_iterator for getting iterator to 
interior/exterior cell +  face index:
std::pair::cell_iterator,  unstigned int>  
get_face_iterator(const unsigned int face_batch_number,
  const unsigned int vector_number, const bool iterior=true,
 const unsigned int fe_component=0)
2) getter for cell_level_index so that the data from get_face_info could be 
used easily (just in case)
 

> That said, for the access of rho and mu on an element, we provide the 
> following access function:
>
>
> https://www.dealii.org/current/doxygen/deal.II/classFEEvaluationBase.html#a442f1ee55e6e80a58e19d53c7184c1a7
>
> This FEEvaluationBase::read_cell_data function expects an 
> AlignedVector> as the argument with as many 
> components as there are cell batches (plus ghosted ones for the MPI case). 
> So in your code, you would fill the rho and mu fields by a cell loop (for 
> MPI, please run it for "n_cells = this->data->n_macro_cells() + 
> this->data->n_ghost_cell_batches()". Then you can access that data both on 
> cells and faces, i.e., with FEEvaluation and FEFaceEvaluation.
>
> The implementation resolves the `interior_cells` and `exterior_cells` 
> fields for you.
>
> Retrieving the `user_index` is a bit more involved as we do not provide 
> access to a face iterator.
>
> You could do it approximately as follows (sitting on a face like in your 
> code below):
>
> const unsigned int interior_cell = face2cell_info.cells_interior[element];
> const Triangulation::cell_iterator my_cell =
>   this->data->get_cell_iterator(interior_cell / 
> VectorizedArray::n_array_elements,
>  interior_cell % 
> VectorizedArray::n_array_elements);
> const Triangulation::face_iterator face_it = 
> my_cell->face(face2cell_info.interior_face_no);
>
> Can you try if that works? If yes, we should add some access functions to 
> MatrixFree. Let us know if you are interested in helping with a pull 
> request.
>
 
Ok, it works (at least on active level). I tested that by comparing 
location of quadrature points (code below). But for implementing 
get_face_iterator I  think that it would be safer to do it using  
cell_level_index (eg. numbering of cell may change )
FEFaceEvaluation 
phi_inner(*(this->data),
  
 true);
FEFaceValues fe_face_values(this->data->get_dof_handler().get_fe(),
this->data->get_face_quadrature(),
UpdateFlags::update_quadrature_points  );

  const unsigned int n_faces=this->data->n_inner_face_batches();
  mu_inner_faces.resize(n_faces);
  for (unsigned int face=0; face::n_array_elements > 
face2cell_info=this->data->get_face_info(face);

  for(unsigned int element=0;
  elementdata->n_active_entries_per_face_batch(face) ;
  ++element){
  const unsigned int interior_cell = face2cell_info.cells_interior[element];
// const unsigned int exterior_cell = 
face2cell_info.cells_exterior[element];

  const typename Triangulation::cell_iterator my_cell =
this->data->get_cell_iterator(interior_cell / 
VectorizedArray::n_array_elements,
   interior_cell % 
VectorizedArray::n_array_elements);
  const typename Triangulation::face_iterator face_it = 
my_cell->face(face2cell_info.interior_face_no);

  fe_face_values.reinit(my_cell, face2cell_info.interior_face_no);
  std::vector > iterator_points = 
fe_face_values.get_quadrature_points ();

  for(unsigned int point_number=1; point_number > point_face 
=phi_inner.quadrature_point (point_number);
double difference =0;
for(unsigned int d=0; d1e-8)
std::cout<<"Diffrence:  " 

[deal.II] MGTransferBlockMatrixFree< dim, Number > ::clear rrmoves number of blocks

[Michał Wichrowski]

Dear all,
I found out that calling
 MGTransferBlockMatrixFree< dim, Number > ::clear
sets the number of blocks to 0/ It is followed directly by build(dofs); 
will thow an exeption: 
Additional information: 
Dimension 0 not equal to 2.
There is no way to set number of block properly other than using 
MGConstrainedDoFs 
 . 
In case of DG where  there are no constrains this requires creating dummy 
objects.
Moreover the default contructor (without parameters) does not work (error 
in compilation - deleted function). Using it could be useful for DG.  I 
would suggest:
1) adding constuctor with number of blocks
2) adding function that sets number of blocks

Best 
Michal

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e7d391cb-5014-4dc7-83a7-b595ac7233ce%40googlegroups.com.

Re: [deal.II] MGTransferBlockMatrixFree< dim, Number > ::clear rrmoves number of blocks

[Michał Wichrowski]

>
>
> So you want to a separate function that calls the current interface with 
> default constructed MGConstraindedDoFs objects?
> We are happy to have a look and discuss if you create a pull request for 
> this.
>
> Thanks, for reporting the problem with the default constructor. I will 
> have a look.
>
> Best,
> Daniel
>
> Ok, I'll try to do it (little later, I'm too busy now).

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1e47bb7e-9a89-4fc3-a261-abbba7b8ffc2%40googlegroups.com.

[deal.II] Adding random filling functionality to vector classes

[Michał Wichrowski]

Dear all, 
It nothing big, but I thing adding method fill_random(const double& a=-1., 
const double &b=1.) to vector classes that will fill a vector with random 
values (from range [a,b]) would be useful. I'm using it in most of my 
programs,  it's great for testing solvers, especially multigrid smoothers. 
Best,
Michal

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/339524fd-6b47-4f2c-8197-8547d4af313d%40googlegroups.com.

[deal.II] Extracting DoF indice from FEEvaluation / Distributing from FEEvaluation to row of a matrix

[Michał Wichrowski]

Dear all,
I need to construct a sparse matrix from MatrixFree operator. I want to do 
it in a similar way to present in step-37 for computing diagonal. The only 
difference will be that the computed values from  FEEValuation and 
FEFaceEvaluation have to be distributed to the specified row of a matrix. 
The code will look like, I need to get local_dof_indices somehow:

// Loop is executed without threads
for (unsigned int cell=0; cell 
(), j);
phi.submit_dof_value(make_vectorized_array(1.), i);
// local integration operator:
integrate_cell(phi,cell);

phi.integrate (false, true);
// Here I have phi filled with a row of the matrix 
// now I need to know global indices of dofs...
// And I have no clue how to obtain local_dof_indices needed below.
for (unsigned int j=0; jdata->n_components_filled(cell); ++v)
matrix(local_dof_indices[i][v],local_dof_indices[j][v] )+=
phi.get_dof_value(j)[v];
}
} 
}

The code for FEFaceVAlues will be very similar.

Moreover, I'm working with a vector-valued problem so local_dof_indices 
will need 3 indices ( local dof number, index in vectorization and 
component of the vector)
It would be great it FEEvaluation had a function for distributing to a 
matrix so that two loops are replaced with phi.distribute_local_to_global 
(matrix, i);
I only need to do that somehow, it does not have to be very efficient. 
There are three main purposes for having that function:
1) coarse matrix for direct solver
2) initialization of preconditioner based on sparse matrix 
3) system matrix for validating code 

Having this functionality would allow me to write a class for an arbitrary 
matrix-free operator that form local integration formula will be able to 
automatically:
1) compute matrix-vector product
2) compute diagonal 
3) compute sparse matrix of operator.

The operator will be specyfied by defining a function :
void integrate_cell(FEEValuation & phi, unsigned int cell);
And similar ones for face and boundary in case of DG. I think it will be 
possible for both vector and scalar based 2nd order operators. Also 
non-square operator could be possible. And using this, the block operators 
may be constructed (eg. Stokes).

Best,
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e3f68a4b-4361-4167-9514-d89e746bfd8b%40googlegroups.com.

[deal.II] Bug:Chebyshev preconditioner of 0th degree

[Michał Wichrowski]

Dear all,
The documentation of Chebyshev precondition says that:

"Note that if the degree variable is set to *zero*, the Chebyshev iteration 
corresponds to a Jacobi preconditioner"

while actually setting degree to 0 triggers Assertion at line 2211 of 
precondition.h. Moreover, the PreconditionChebyshev works well with 
degree=0 in release mode, so I think removing assetion will resolve the 
problem.

If Chebyshev preconditioner is not intended to be used with 0th degree this 
should be added to documentation (and documentation in current form is 
missleading)

Best,
MIchaĺ

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/052f0208-e7b2-41aa-8876-c0df06ca4a92%40googlegroups.com.

Re: [deal.II] Bug:Chebyshev preconditioner of 0th degree

[Michał Wichrowski]

W dniu sobota, 12 października 2019 14:29:23 UTC+2 użytkownik Daniel Arndt 
napisał:
>
> MIchaĺ,
>
> The documentation in the developer version or the current version (
> https://www.dealii.org/developer/doxygen/deal.II/classPreconditionChebyshev.html)
>  
> says:
> "[...] Note that if the degree variable is set to one, the Chebyshev 
> iteration corresponds to a Jacobi preconditioner (or the underlying 
> preconditioner type) with relaxation parameter according to the specified 
> smoothing range. [...]"
> and the assertion is at line 2211 only in the developer version. The text 
> you are referring to can be found in version 9.0.0.
>
> Best,
> Daniel
>
>>
>> You're right, I was following wrong version of documentation. 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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/576b0ce7-5eb4-4d9a-af2a-80262ce3d50e%40googlegroups.com.

[deal.II] Matrix-free DG with adaptive refinements

[Michał Wichrowski]

Dear all,
I modified step-59 function run() so that in only makes grid and does local 
refinement and initializes MatrixFree object.  This triggers an assertion 
when calling
mg_mf_storage_level->reinit(dof_handler,
dummy,
QGauss<1>(fe.degree + 1),
additional_data);
in system_setup().

Function run():
Point upper_right;
upper_right[0] = 2.5;
for (unsigned int d = 1; d < dim; ++d)
  upper_right[d] = 2.8;
GridGenerator::hyper_rectangle(triangulation,
   Point(),
   upper_right);
triangulation.begin_active()->face(0)->set_boundary_id(10);
triangulation.begin_active()->face(1)->set_boundary_id(11);
triangulation.begin_active()->face(2)->set_boundary_id(0);
for (unsigned int f = 3; f < GeometryInfo::faces_per_cell; 
++f)
  triangulation.begin_active()->face(f)->set_boundary_id(1);

std::vector::cell_iterator>>
  periodic_faces;
GridTools::collect_periodic_faces(
  triangulation, 10, 11, 0, periodic_faces);
triangulation.add_periodicity(periodic_faces);

triangulation.refine_global(1);

Vector 
estimated_error_per_cell(triangulation.n_active_cells());
for (unsigned int i = 0; i < estimated_error_per_cell.size(); 
++i)
  estimated_error_per_cell[i] = std::sin(i) * std::sin(i);


parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
triangulation, estimated_error_per_cell, 0.3, 0.0);
triangulation.execute_coarsening_and_refinement();
  
setup_system();
  

Error message:


An error occurred in line <4993> of file 
 in function
void dealii::FESubfaceValues::reinit(const typename 
dealii::Triangulation::cell_iterator&, unsigned int, 
unsigned int) [with int dim = 3; int spacedim = 3; typename 
dealii::Triangulation::cell_iterator = 
dealii::TriaIterator >]
The violated condition was: 
subface_no < cell->face(face_no)->n_children()
Additional information: 
Index 0 is not in the half-open range [0,0). In the current case, this 
half-open range is in fact empty, suggesting that you are accessing an 
element of an empty collection such as a vector that has not been set to 
the correct size.

Stacktrace:
---

I have no idea what may be causing this. 

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/69586e52-7e99-41e4-a7f1-930cffd6c046%40googlegroups.com.

[deal.II] Matrix-free: Continuous elements and boundary (face) integrals

[Michał Wichrowski]

Deal all,
I want to weakly impose boundary conditions on, let's say, Laplace problem 
(in fact I'm working on Stokes).  I'm considering [FE_Q<2>(2)^2]   ( that 
is   FESystem<2>[FE_Q<2>(2)^2] + FE_Q<2>(1) for Stokes).
As a MatrixFree loop I'm using:
  this->data->loop(&SymGradMass::local_apply, &SymGradMass::apply_face,
   &SymGradMass::apply_boundary,
   this, dst, src,
   /*zero_dst =*/false,
   MatrixFree::DataAccessOnFaces::gradients,
   MatrixFree::DataAccessOnFaces::gradients);

where apply_face function is empty (only local_apply and apply_boundary 
perform integration ).

 After initializing MatrixFree object with following data:

typename MatrixFree::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree::AdditionalData::none;
additional_data.mapping_update_flags =
  (update_gradients | update_JxW_values | update_quadrature_points);
additional_data.mapping_update_flags_boundary_faces =
(update_gradients | update_JxW_values | update_normal_vectors |
 update_quadrature_points);
additional_data.mapping_update_flags_inner_faces =
(update_default);

I get end error (at the end of the email).  It looks like the MatrixFree is 
trying to compute coupling of DoF over face but on the other side that 
information is not provided (by FE?). I suppose this only requires simple 
fix:  in case if 
mapping_update_flags_inner_faces =update_default
proceed as in standard case (no boundary integrals) - the boundary 
integrals does not generate additional DoF coupling.

Best,
Michał Wichrowski

An error occurred in line <786> of file 
 in function
dealii::types::global_dof_index 
dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const
The violated condition was: 
static_cast::type)>:: type>(local_index) < static_cast::type)>:: type>(local_size() + n_ghost_indices_data)
Additional information: 
Index 4294967295 is not in the half-open range [0,12832).

Stacktrace:
---
#0  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const
#1  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void 
dealii::MatrixFree<2, double, dealii::VectorizedArray 
>::initialize_indices(std::vector 
const*, std::allocator const*> > const&, 
std::vector > const&, 
dealii::MatrixFree<2, double, dealii::VectorizedArray 
>::AdditionalData const&)
#2  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void 
dealii::MatrixFree<2, double, dealii::VectorizedArray 
>::internal_reinit(dealii::Mapping<2, 2> const&, 
std::vector const*, 
std::allocator const*> > const&, 
std::vector const*, 
std::allocator const*> > const&, 
std::vector > const&, 
std::vector, 
std::allocator > > const&, dealii::MatrixFree<2, 
double, dealii::VectorizedArray >::AdditionalData const&)
#3  ./fsi_matrix_free_2d.g: void dealii::MatrixFree<2, double, 
dealii::VectorizedArray >::reinit, 
dealii::QGauss<1>, double>(dealii::Mapping<2, 2> const&, 
std::vector const*, 
std::allocator const*> > const&, 
std::vector const*, 
std::allocator const*> > const&, 
dealii::QGauss<1> const&, dealii::MatrixFree<2, double, 
dealii::VectorizedArray >::AdditionalData const&)
#4  ./fsi_matrix_free_2d.g: StokesSolver<2, 2>::setup_system()
#5  ./fsi_matrix_free_2d.g: StokesSolver<2, 
2>::initialize(std::shared_ptr >&)
#6  ./fsi_matrix_free_2d.g: FSIInterface<2, 
2>::FSIInterface(RunTimeParameters const&)
#7  ./fsi_matrix_free_2d.g: main



-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5cabfdc7-7db3-484d-b705-e31be5c4bd68%40googlegroups.com.

[deal.II] Assembling sparse matrix from Matrix-free vmult with constrains

[Michał Wichrowski]

Dear  all,
I've got matrix-free multigrid solver for Stokes problem. The main 
bottleneck is solution of coarse problem, so I tried to assemble the 
regular sparse matrix and use direct solver. Since the coarse problem is 
(relatively) small, I used vmults by unit vector to obtain columns of the 
matrix.  This is my code:


setup_matrix(); // generates sparsity patters and reinits matrices

LinearAlgebra::distributed::BlockVector  dst;
LinearAlgebra::distributed::BlockVector  src;

src.reinit(2);
for (unsigned int b = 0; b < 2; ++b)
stokes_matrix.get_matrix_free()->initialize_dof_vector(
src.block(b), b);
src.collect_sizes();
src =0;


dst.reinit(2);
dst.block(0).reinit(owned_dofs_u, relevant_dofs_u, MPI_COMM_WORLD);
dst.block(1).reinit(owned_dofs_p, relevant_dofs_p, MPI_COMM_WORLD);
dst.collect_sizes();

for(types::global_dof_index i =0; i< owned_dofs_u.size(); ++i){
src=0;
dst=0;
if(owned_dofs_u.is_element(i) )
src.block(0)(i)=1;

src.compress(VectorOperation::insert);
stokes_matrix.vmult(dst, src);

dst.update_ghost_values();

for(IndexSet::ElementIterator  index_iter =owned_dofs_u.begin();
index_iter != owned_dofs_u.end();
++index_iter){
if(dst.block(0)(*index_iter)!=0 )
block_A->set(*index_iter,i, dst.block(0)(*index_iter) );
}

}

block_A->compress(VectorOperation::unknown);

Without constrains the matrix matches the matrix-free operator, but with 
constrains present it does not. What is the proper way to assemble the 
matrix with vmult?

Best,
Michał Wichrowski



-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c1bf23f7-a8fb-41cb-b24c-877afc1b16b1%40googlegroups.com.

[deal.II] Re: Assembling sparse matrix from Matrix-free vmult with constrains

[Michał Wichrowski]

It cames out it is not so easy to provide simplest working example. I've 
tried to remake step-37, but surprisingly I was not able to reproduce the 
problem. I'll take a closer look what is causing a problem. 

W dniu wtorek, 10 marca 2020 07:46:11 UTC+1 użytkownik peterrum napisał:
>
> Hi Michal,
>
> any chance that you post or send me a small runnable test program.
>
> By the way, there is an open PR (
> https://github.com/dealii/dealii/pull/9343) for the computation of the 
> diagonal in a matrix-free manner. Once this is merged, I will work the 
> matrix-free assembly of sparse matrices.
>
> Thanks,
> Peter
>
> On Monday, 9 March 2020 15:26:10 UTC+1, Michał Wichrowski wrote:
>>
>> Dear  all,
>> I've got matrix-free multigrid solver for Stokes problem. The main 
>> bottleneck is solution of coarse problem, so I tried to assemble the 
>> regular sparse matrix and use direct solver. Since the coarse problem is 
>> (relatively) small, I used vmults by unit vector to obtain columns of the 
>> matrix.  This is my code:
>>
>>
>> setup_matrix(); // generates sparsity patters and reinits matrices
>>
>> LinearAlgebra::distributed::BlockVector  dst;
>> LinearAlgebra::distributed::BlockVector  src;
>>
>> src.reinit(2);
>> for (unsigned int b = 0; b < 2; ++b)
>> stokes_matrix.get_matrix_free()->initialize_dof_vector(
>> src.block(b), b);
>> src.collect_sizes();
>> src =0;
>>
>>
>> dst.reinit(2);
>> dst.block(0).reinit(owned_dofs_u, relevant_dofs_u, MPI_COMM_WORLD);
>> dst.block(1).reinit(owned_dofs_p, relevant_dofs_p, MPI_COMM_WORLD);
>> dst.collect_sizes();
>>
>> for(types::global_dof_index i =0; i< owned_dofs_u.size(); ++i){
>> src=0;
>> dst=0;
>> if(owned_dofs_u.is_element(i) )
>> src.block(0)(i)=1;
>>
>> src.compress(VectorOperation::insert);
>> stokes_matrix.vmult(dst, src);
>>
>> dst.update_ghost_values();
>>
>> for(IndexSet::ElementIterator  index_iter =owned_dofs_u.begin();
>> index_iter != owned_dofs_u.end();
>> ++index_iter){
>> if(dst.block(0)(*index_iter)!=0 )
>> block_A->set(*index_iter,i, dst.block(0)(*index_iter) );
>> }
>>
>> }
>>
>> block_A->compress(VectorOperation::unknown);
>>
>> Without constrains the matrix matches the matrix-free operator, but with 
>> constrains present it does not. What is the proper way to assemble the 
>> matrix with vmult?
>>
>> Best,
>> Michał Wichrowski
>>
>>
>>
>>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/3944b062-fc59-4894-bc89-c43e017c16e9%40googlegroups.com.

Re: [deal.II] Assembling sparse matrix from Matrix-free vmult with constrains

[Michał Wichrowski]

W dniu wtorek, 10 marca 2020 17:33:43 UTC+1 użytkownik Wolfgang Bangerth 
napisał:
>
> On 3/9/20 8:26 AM, Michał Wichrowski wrote: 
> > I've got matrix-free multigrid solver for Stokes problem. The main 
> > bottleneck is solution of coarse problem, so I tried to assemble the 
> > regular sparse matrix and use direct solver. Since the coarse problem is 
> > (relatively) small, I used vmults by unit vector to obtain columns of 
> > the matrix. 
>
> This is unrelated to your original question, but worth mentioning 
> anyway: Conceptually, you are multiplying a unit vector by a sparse 
> matrix, and what you should get is a vector with nonzero entries in only 
> those rows for which the matrix has a nonzero entry in the column 
> corresponding to the unit vector. 
>
> This is inefficient, because you will have to do N such matrix-vector 
> products. But if you know which entries of the matrix are non-zero 
> (which you do, because you know which DoFs couple with each other), then 
> you can multiply with a sum of unit vectors, get an output vector, and 
> you will know which entries of the output vector correspond to the 
> product with each of these unit vectors. You should be able to get the 
> number of matrix-vector products down from N to N/100 or N/1000 for 
> sufficiently large problems -- i.e., *much much faster*. 
>
> (In fact, I think that we should at one point add a function that takes 
> an operator and that computes a matrix representation. If we would also 
> give that function a SparsityPattern, it could do the algorithm outlined 
> above automatically. As always, any help would be welcome!) 
>

Yes, I thought about that, but I my case is not worth an effort. I'm doing 
it to assemble coarse matrix, so the problem is relatively small (N_dof 
around 500). The assemble time is  unnoticeable comparing to solution time, 
even in case of small problems. Moreover, I need to apply inverse of 
some-kind Schur complement matrix (on coarse level only), and for that it 
would be extremely hard to guess sparsity pattern (it is still sparse). 
Fortunately, this matrix is even smaller. 

Your strategy (if implemented correctly) would require constant number of 
matrix-vector multiplications so the complexity will be linear. The number 
of vmults will depend on the mesh connectivity and type of FE. So the 
overall complexity will  be linear with the problem size (instead (O(N^2), 
as I understand from Your email). But indeed, it will be much faster. 
Best,
Michał

>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: 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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6992d3bb-ad0e-4c96-b539-f2b5db84cf21%40googlegroups.com.

[deal.II] Re: Matrix-free: Continuous elements and boundary (face) integrals

[Michał Wichrowski]

I tried to fix that by myself, but I could not locate the source code of 
MatrixFree::internal_reinit

W dniu poniedziałek, 9 marca 2020 14:20:16 UTC+1 użytkownik Michał 
Wichrowski napisał:
>
> Deal all,
> I want to weakly impose boundary conditions on, let's say, Laplace problem 
> (in fact I'm working on Stokes).  I'm considering [FE_Q<2>(2)^2]   ( that 
> is   FESystem<2>[FE_Q<2>(2)^2] + FE_Q<2>(1) for Stokes).
> As a MatrixFree loop I'm using:
>   this->data->loop(&SymGradMass::local_apply, &SymGradMass::apply_face,
>&SymGradMass::apply_boundary,
>this, dst, src,
>/*zero_dst =*/false,
>MatrixFree::DataAccessOnFaces::gradients,
>MatrixFree::DataAccessOnFaces::gradients);
>
> where apply_face function is empty (only local_apply and apply_boundary 
> perform integration ).
>
>  After initializing MatrixFree object with following data:
>
> typename MatrixFree::AdditionalData additional_data;
> additional_data.tasks_parallel_scheme =
> MatrixFree::AdditionalData::none;
> additional_data.mapping_update_flags =
>   (update_gradients | update_JxW_values | update_quadrature_points);
> additional_data.mapping_update_flags_boundary_faces =
> (update_gradients | update_JxW_values | update_normal_vectors |
>  update_quadrature_points);
> additional_data.mapping_update_flags_inner_faces =
> (update_default);
>
> I get end error (at the end of the email).  It looks like the MatrixFree 
> is trying to compute coupling of DoF over face but on the other side that 
> information is not provided (by FE?). I suppose this only requires simple 
> fix:  in case if 
> mapping_update_flags_inner_faces =update_default
> proceed as in standard case (no boundary integrals) - the boundary 
> integrals does not generate additional DoF coupling.
>
> Best,
> Michał Wichrowski
>
> An error occurred in line <786> of file 
>  in function
> dealii::types::global_dof_index 
> dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const
> The violated condition was: 
> static_cast std::common_type n_ghost_indices_data)>::type)>:: type>(local_index) < static_cast ::dealii::internal::argument_type std::common_type n_ghost_indices_data)>::type)>:: type>(local_size() + n_ghost_indices_data)
> Additional information: 
> Index 4294967295 is not in the half-open range [0,12832).
>
> Stacktrace:
> ---
> #0  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const
> #1  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void 
> dealii::MatrixFree<2, double, dealii::VectorizedArray 
> >::initialize_indices(std::vector 
> const*, std::allocator const*> > const&, 
> std::vector > const&, 
> dealii::MatrixFree<2, double, dealii::VectorizedArray 
> >::AdditionalData const&)
> #2  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void 
> dealii::MatrixFree<2, double, dealii::VectorizedArray 
> >::internal_reinit(dealii::Mapping<2, 2> const&, 
> std::vector const*, 
> std::allocator const*> > const&, 
> std::vector const*, 
> std::allocator const*> > const&, 
> std::vector > const&, 
> std::vector, 
> std::allocator > > const&, dealii::MatrixFree<2, 
> double, dealii::VectorizedArray >::AdditionalData const&)
> #3  ./fsi_matrix_free_2d.g: void dealii::MatrixFree<2, double, 
> dealii::VectorizedArray >::reinit, 
> dealii::QGauss<1>, double>(dealii::Mapping<2, 2> const&, 
> std::vector const*, 
> std::allocator const*> > const&, 
> std::vector const*, 
> std::allocator const*> > const&, 
> dealii::QGauss<1> const&, dealii::MatrixFree<2, double, 
> dealii::VectorizedArray >::AdditionalData const&)
> #4  ./fsi_matrix_free_2d.g: StokesSolver<2, 2>::setup_system()
> #5  ./fsi_matrix_free_2d.g: StokesSolver<2, 
> 2>::initialize(std::shared_ptr >&)
> #6  ./fsi_matrix_free_2d.g: FSIInterface<2, 
> 2>::FSIInterface(RunTimeParameters const&)
> #7  ./fsi_matrix_free_2d.g: main
> 
>
>
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/469c66f0-f9b1-4bc4-955f-8f6b5e19662c%40googlegroups.com.

Re: [deal.II] Assembling sparse matrix from Matrix-free vmult with constrains

[Michał Wichrowski]

It looks like it was hidden bug inside my code. I've wrote the coarse grid 
direct solver for matrix-free operators. It ignores sparsity pattern and 
assembles the matrix by doing N vmults. The interface may be extended for 
optimalization. I think it might be an useful tool. Usage, for step-37:
constructor needs vector partitioning:
  CoarseDirectSolver coarse_solver(owned_dofs_u,relevant_dofs_u );
This  
  coarse_solver.initialize(mg_matrices[0]/*, dof_handler*/);
or
  coarse_solver.initialize(mg_matrices[0], dof_handler);

assembles the marix. In second case the matrix is initilized with sparsity 
pattern, it could be optimized to reduce the number of vmults.

Michał

W dniu wtorek, 10 marca 2020 18:28:15 UTC+1 użytkownik Michał Wichrowski 
napisał:
>
>
>
> W dniu wtorek, 10 marca 2020 17:33:43 UTC+1 użytkownik Wolfgang Bangerth 
> napisał:
>>
>> On 3/9/20 8:26 AM, Michał Wichrowski wrote: 
>> > I've got matrix-free multigrid solver for Stokes problem. The main 
>> > bottleneck is solution of coarse problem, so I tried to assemble the 
>> > regular sparse matrix and use direct solver. Since the coarse problem 
>> is 
>> > (relatively) small, I used vmults by unit vector to obtain columns of 
>> > the matrix. 
>>
>> This is unrelated to your original question, but worth mentioning 
>> anyway: Conceptually, you are multiplying a unit vector by a sparse 
>> matrix, and what you should get is a vector with nonzero entries in only 
>> those rows for which the matrix has a nonzero entry in the column 
>> corresponding to the unit vector. 
>>
>> This is inefficient, because you will have to do N such matrix-vector 
>> products. But if you know which entries of the matrix are non-zero 
>> (which you do, because you know which DoFs couple with each other), then 
>> you can multiply with a sum of unit vectors, get an output vector, and 
>> you will know which entries of the output vector correspond to the 
>> product with each of these unit vectors. You should be able to get the 
>> number of matrix-vector products down from N to N/100 or N/1000 for 
>> sufficiently large problems -- i.e., *much much faster*. 
>>
>> (In fact, I think that we should at one point add a function that takes 
>> an operator and that computes a matrix representation. If we would also 
>> give that function a SparsityPattern, it could do the algorithm outlined 
>> above automatically. As always, any help would be welcome!) 
>>
>
> Yes, I thought about that, but I my case is not worth an effort. I'm doing 
> it to assemble coarse matrix, so the problem is relatively small (N_dof 
> around 500). The assemble time is  unnoticeable comparing to solution time, 
> even in case of small problems. Moreover, I need to apply inverse of 
> some-kind Schur complement matrix (on coarse level only), and for that it 
> would be extremely hard to guess sparsity pattern (it is still sparse). 
> Fortunately, this matrix is even smaller. 
>
> Your strategy (if implemented correctly) would require constant number of 
> matrix-vector multiplications so the complexity will be linear. The number 
> of vmults will depend on the mesh connectivity and type of FE. So the 
> overall complexity will  be linear with the problem size (instead (O(N^2), 
> as I understand from Your email). But indeed, it will be much faster. 
> Best,
> Michał
>
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth  email: 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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f8204afb-1750-4459-a7d1-67c33e4a256e%40googlegroups.com.
/*
 * assemble_coarse_matrix.h
 *
 *  Created on: Mar 10, 2020
 *  Author: mwichro
 */

#ifndef ASSEMBLE_COARSE_MATRIX_H_
#define ASSEMBLE_COARSE_MATRIX_H_


#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 


#include 
#include 
#include 
#include 
#include 



using namespace dealii;

[deal.II] Coarse direct solver for MatrixFree (also block version)

[Michał Wichrowski]

Dear all,
I've wrote a interface for Trilinos direct solver so that it may be used as 
coarse solver for MG. It assembles the matrix by using N vmults (so this 
part is very inefficient), so it will only pay off if solver is used 
multiple times.  I've got also block version. It is possible to extend 
block to work as an (efficient) direct solver for block matrices. Since we 
do not have a parallel direct solver for block matrices it might be an 
useful tool.  

Best,
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f28fea9d-64ab-4422-a2b6-66bc28aab5d4%40googlegroups.com.
/*
 * block_direct_solver.h
 *
 *  Created on: Mar 18, 2020
 *  Author: mwichro
 */

#ifndef INCLUDE_BLOCK_DIRECT_SOLVER_H_
#define INCLUDE_BLOCK_DIRECT_SOLVER_H_

#include 
#include 

#include 

#include 
#include 
#include 
#include 

#include 


#include 

using namespace dealii;

template 
class BlockDirectSolver
:
	public MGCoarseGridBase>{

public:
friend class BlockCoarseTest;
	typedef LevelNumber number;
	typedef LinearAlgebra::distributed::BlockVector BlockVectorType;


	//Constructor. Requires owned_partitioning[block][MPI process],
	// relevant_partitioning[block] and MPI communicator.
	// The owned partitioning for each block may be obtained by using
	// DoFHandler::compute_locally_owned_dofs_per_processor() or in case of MG
	// DoFHandler::compute_locally_owned_mg_dofs_per_processor(level)


	BlockDirectSolver(const std::vector>& owned_partitioning_,
			const std::vector& relevant_dofs_,
			const MPI_Comm & communicator);

	//Computes TrilinosWrappers::SparseMatrix using internal DoF renumbering,
	// initializes direct solver (TrilinosWrappers::SolverDirect)
	// Level will be checked if  operator() is used
	template
	void initialize(const Operator & vmult_operator,
			const unsigned int& lvl=0);
//
//	template
//		void initialize(const Operator & vmult_operator,
//const std::vector*> &dof_handlers,
//const unsigned int& lvl=0);
//


	void vmult(BlockVectorType &dst,
			const BlockVectorType &src) const;
	//Interface to act as a coarse grid solver.
	//Applies inverse of operator used in initialize.
	virtual void
	  operator()(const unsigned int lvl,
			  BlockVectorType &  dst,
	 const BlockVectorType &src) const;
private:


	typedef TrilinosWrappers::MPI::Vector InternalVectorType;

	const MPI_Comm  mpi_comm;
	const unsigned int this_mpi_process;
	const unsigned int n_mpi_process;

	int level;
	const unsigned int n_blocks;

	types::global_dof_index n_local_dofs;

	std::vector owned_dofs;
	//[block][process]


	const std::vector> owned_partitioning;
	const std::vector relevant_dofs;

	//Searching for owning processor is in order specified by following list
	mutable std::list searching_order;

	//Shitfs of dof indices, used for internal renumbering
	std::vector> dof_shifts;

	std::vector first_block_index;
	std::vector last_block_index;

	//IndexSet of reordered owned dofs.
	IndexSet shifted_owned;

	std::unique_ptr  matrix;
	SolverControl solver_control;
	mutable TrilinosWrappers::SolverDirect inverse;


	//Computes internal DoF number from external (dof, block) numbering.
	//First checks if dof is locally owned, then searches other index sets, in order
	// specified by searching_order.
	types::global_dof_index ext2int(const types::global_dof_index& external_dof,
			const unsigned int block) const ;


	//Not implemented, not needed. Placeholder.
	 std::pair  int2ext (types::global_dof_index &internal_dof) const;

};

template 
BlockDirectSolver::BlockDirectSolver(const std::vector>& owned_partitioning_,
		const std::vector& relevant_dofs_,
		const MPI_Comm & communicator)
		:
		mpi_comm(communicator)
	 ,  this_mpi_process(Utilities::MPI::this_mpi_process(communicator))
	 ,  n_mpi_process(Utilities::MPI::n_mpi_processes(mpi_comm))
	 ,	level(-1)
	 ,  n_blocks(relevant_dofs_.size())
	 ,  n_local_dofs(0)
	 ,  owned_partitioning(owned_partitioning_)
	 ,  relevant_dofs(relevant_dofs_)
	 ,  solver_control(10, 1e-14)
	 ,  inverse(solver_control)
{

	AssertDimension(owned_partitioning.size(), n_blocks);
	owned_dofs.resize(n_blocks);
	for(unsigned int b=0; b starts(n_mpi_process,0);
	types::global_dof_index n_dofs=0;




	for(unsigned int i=1; i< n_mpi_process; ++i){
		starts[i]=starts[i-1];

		for(unsigned int b=0; b <  n_blocks; ++b){
			starts[i]+=owned_partitioning[b][i-1].n_elements();

		}
	}


	for(unsigned int b=0; b <  n_blocks; ++b)
			n_dofs+=owned_partitioning[b][this_mpi_process].size();


	for(unsigned int p=0; p
template
void
BlockDirectSolver::initialize(const Operator & vmult_operator

[deal.II] MtrixFree for biharmonic problem: submit_laplacian

[Michał Wichrowski]

Dear all,

I was trying to implement MatrixFree biharmonic solver. We already have 
most of functionality needed to do it (DG,  computing laplacian on both 
faces and cells), but the one thing missing is extending FEEvaluation and 
FEFaceEvaluation to support integration of 2nd order formulae: 
submit_laplacian() and extending integrate() would resolve the problem and 
allow matrix-free version  of step-47. 

The biharmonic solver may be also useful in other parts of library, 
especially handling mesh deformation and transforms.

Best,
Michał
 

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2448bf19-ee2f-4bdb-8ddd-f192d064db0e%40googlegroups.com.

Re: [deal.II] MtrixFree for biharmonic problem: submit_laplacian

[Michał Wichrowski]

I've just looked through the possible solution strategies and what could be 
done inside existing deal.II functionality. While implementing I noticed 
that there is missing option of subimiting laplacian and I thought that 
will be something easy to do. Now I see where the problem is.

I need a quick fix for my FSI solver (moving mesh) so I'm rewriting step-47 
in parallel with multigrid solver. It required some fixes in deal.ii but 
think it is going to work. I'll post the patch when it will be working. 
Thanks,
Michał

W dniu sobota, 25 kwietnia 2020 14:47:34 UTC+2 użytkownik Martin 
Kronbichler napisał:
>
> Dear Michal,
>
> That is great news. There is interest in this functionality also by Julius 
> Witte (I added him here) and there were some steps done.
>
> The main difficulty in terms of user interface is the fact that the 
> submit_laplacian() is a bit nastier than the other functions: on non-affine 
> geometries, the Laplacian (or Hessian) of the test function in real space 
> gives contributions to both the Laplacian and the gradient in reference 
> space. So we either need a `submit_laplacian()` function that both inserts 
> the desired contribution to test by the gradient and the Laplacian - which 
> breaks symmetry with the get_gradient()/get_laplacian() methods - or we 
> need to have a second field that stores the accumulated part for the 
> reference-cell gradient of the test function and merges this inside 
> 'integrate'.
>
> When you say "we already have most functionality needed to to it (DG, 
> computing laplacian on both faces and cells)", do you imply that you have 
> already started with the work needed for get_laplacian() on faces? At least 
> on dealii master, we lack both the ability to evaluate the second 
> derivatives in tangential directions in reference space and the Jacobian 
> gradient needed for non-affine cell shapes.
>
> In any case, I would be happy to help and join the efforts, and I think so 
> is Julius. We have an imminent release in a few weeks for deal.II 
> (scheduled feature freeze May 8), so the goal would be to try to achieve 
> this soon after the release.
>
> Best,
> Martin
> On 25.04.20 12:28, Michał Wichrowski wrote:
>
> Dear all,
>
> I was trying to implement MatrixFree biharmonic solver. We already have 
> most of functionality needed to do it (DG,  computing laplacian on both 
> faces and cells), but the one thing missing is extending FEEvaluation and 
> FEFaceEvaluation to support integration of 2nd order formulae: 
> submit_laplacian() and extending integrate() would resolve the problem and 
> allow matrix-free version  of step-47. 
>
> The biharmonic solver may be also useful in other parts of library, 
> especially handling mesh deformation and transforms.
>
> Best,
> Michał
>  
> -- 
> 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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/2448bf19-ee2f-4bdb-8ddd-f192d064db0e%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/2448bf19-ee2f-4bdb-8ddd-f192d064db0e%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1c83921e-9659-47b7-9a5a-c908e94543a6%40googlegroups.com.

[deal.II] MappingFEField with MatrixFree: MG matrices does not use proper mappling

[Michał Wichrowski]

Dear all,
I've got multigrid solver matrix-free for FSI in ALE formulation. For 
deforming grid I'm using MappingFEField, stored 
in std::shared_ptr>. On the coarsest level I'm assembling the 
sparse matrix in standard way, I then initialize direct solver (MUMPS from 
Trilinos) and than I'm using it as a preconditioner 
in MGCoarseGridIterativeSolver. So I should get convergance of coarse 
solver in 1 iteration (or at most several iterations if the matrix is 
really badly conditioned). I also have a test that check if the matrix free 
operator is the same as assembled matrix (by doing  N vmults).

In case when the mesh is undeformed this is exactly what happens i.e, the 
solver coverges in 1 iteration. However if mesh is deformed the number of 
iterations of coarse level grows. The number of coarse iterations in case 
of Turek FSI3 benchmark may exceed 30. 

I've double-checked the solver, I'm updating mapping on every time step:
  system_mf_storage->update_mapping(*mapping);
  system_mf_storage_no_dirichlet->update_mapping(*mapping);
  for (unsigned int lvl = 0; lvl < triangulation.n_global_levels(); ++lvl)
mg_level_data[lvl].update_mapping(*(level_mappings[lvl])); 
//level_mappings are shared pointers, the same as mapping.
I am also sure that I'm interpolating fine mapping to mg levels. 

The problem occurs only on level matrices, the system matrix that I'm 
solving is OK. The coarse solver assembler gets the same mapping as 
MatrixFree object. 

I've only spotted the problem at coarsest level, but I think it is also an 
issue on other levels. I am using 9.3.0-pre (cluster) and 9.2.0-pre (my 
PC)  versions of Deal.II, the same problem appears on both.

Best ,
Michał

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d9cb2f8a-89bc-4c21-9ce0-19c3085d05f2o%40googlegroups.com.

[deal.II] MGConstrainedDoFs::initialize with no constrains

[Michał Wichrowski]

Dear all,
I'm considering a multigrid method for the Stokes problem. I need to have 
an approximation of Schur complement at each level that is represented by a 
matrix-free operator.

 The problem is that in case of adaptive mesh refinement  MGConstrainedDoFs 
assumes zero boundary condition at refinement edges. In my case that 
behaviour is unwanted. I have looked through the implementation of 
MGConstrainedDoFs and it looks like it will require changing initialize() 
function by adding 1 additional parameter:

MGConstrainedDoFs::initialize ( const DoFHandler< dim, spacedim > & dof,
const MGLevelObject< IndexSet > & level_relevant_dofs = 
MGLevelObject() 
const bool& set_zero_constrain_at_refinement_edge=true) 

And adding if in line 331:
if(set_zero_constrain_at_refinement_edge)
  MGTools::extract_inner_interface_dofs(dof, refinement_edge_indices);

The other option to implement my MG preconditioner is to 
rewrite MGConstrainedDoFs together with all MG classes that rely on it.

Best,
Michał


-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/8c43a6e1-875e-4150-b15d-14b56ab438dan%40googlegroups.com.

[deal.II] Setting Eclipse with CMake

[Michał Wichrowski]

Dear all,
The Cmake generator of Eclipse project does not work very well, it was not 
developed for years. As a consequence, the Eclipse project generated by 
CMake (as recommended in dealii Wiki) sometimes results with annoying 
issues with indexer. 
There is a useful plugin: cmake4eclipse that allows eclipse to index a 
project based on CMakeLists. Here are the instructions how to set a project 
using it:
https://stackoverflow.com/a/38716337
It works very well for me. I also have dealii as a project in eclipse so 
I've specyfie
I think including it in the wiki could be useful. 
Additionaly adding info about support for c++11  in Eclipse would be useful:
https://stackoverflow.com/a/24627391
(for example without it C++ shared/unique pointers are not indexed)

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e537c866-5a6d-40e1-9e78-a60eebc736ccn%40googlegroups.com.

Re: [deal.II] Setting Eclipse with CMake

[Michał Wichrowski]

Yes, some years ago Wolfgang gave me rights to modify wiki, but I didn't 
want to do it without permission.

środa, 3 lutego 2021 o 10:19:22 UTC+1 luca@gmail.com napisał(a):

> Dear Michal, 
>
> can you modify the wiki yourself? It should be possible to open a PR with 
> the suggested modifications. Everyone is welcome to contribute, and I think 
> your information would be very valuable to other users as well.
>
> Best,
> Luca.
>
> > On 30 Jan 2021, at 16:34, Michał Wichrowski  wrote:
> > 
> > Dear all,
> > The Cmake generator of Eclipse project does not work very well, it was 
> not developed for years. As a consequence, the Eclipse project generated by 
> CMake (as recommended in dealii Wiki) sometimes results with annoying 
> issues with indexer. 
> > There is a useful plugin: cmake4eclipse that allows eclipse to index a 
> project based on CMakeLists. Here are the instructions how to set a project 
> using it:
> > https://stackoverflow.com/a/38716337
> > It works very well for me. I also have dealii as a project in eclipse so 
> I've specyfie
> > I think including it in the wiki could be useful. 
> > Additionaly adding info about support for c++11 in Eclipse would be 
> useful:
> > https://stackoverflow.com/a/24627391
> > (for example without it C++ shared/unique pointers are not indexed)
> > 
> > -- 
> > 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/e537c866-5a6d-40e1-9e78-a60eebc736ccn%40googlegroups.com
> .
>
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6b6086c7-2e0c-4b25-81c5-ebaf818accd0n%40googlegroups.com.

Re: [deal.II] How to refine mesh in different field

[Michał Wichrowski]

Jest overwrite error estimate in one domain by zeros:

unsigned int cell_n = 0;
for (const auto & cell : triangulation.active_cell_iterators ){
if (cell-> material_id() == 0 )
estimated_error_per_cell(i) = 0;
++i;
}

sobota, 3 kwietnia 2021 o 16:48:48 UTC+2 hkch...@gmail.com napisał(a):

> Dear deal.II community,
>
> In a fluid-structure interaction simulation, Computational domain was 
> divided into fluid domain and solid domain, for example, cell->material_id() 
> == 0 and cell->material_id() == 1(the way that deal.ii presents which 
> field the cell belong to) respectively, Now I want to refine the mesh in 
> fluid field and that in solid field separately in a fixed time interval.
>
> The tutorial guides, for example, tutorial step-40, show follow example 
> code:
>
> template 
> void LaplaceProblem::refine_grid()
> {
> TimerOutput::Scope 
>  
> t(computing_timer, "refine");
> Vector 
>  
> estimated_error_per_cell(triangulation.n_active_cells 
> 
> ());
> KellyErrorEstimator::estimate 
> 
> (
> dof_handler,
> QGauss 
> (fe.degree 
> + 1),
> std::map , 
> const Function 
>  *>(),
> locally_relevant_solution,
> estimated_error_per_cell);
> parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number 
> 
> (
> triangulation, estimated_error_per_cell, 0.3, 0.03);
> triangulation.execute_coarsening_and_refinement 
> 
> ();
> }
>
> I think the code refines on the whole mesh domain, for my purpose, how to 
> change the code to meet me?
> Could anyone help me out?
>
> Sincerely,
> Chen
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/136c7e2d-2cdc-4abb-b522-965906274eabn%40googlegroups.com.