Re: [deal.II] Re: Set up boundary id in gmsh

[Bruno Blais]

The pleasure is mine.
You will find that people the deal.II community are very helpfull :)!


On Wednesday, 5 August 2020 22:44:11 UTC-4, yuesu jin wrote:
>
> Dear Bruno,
>   Thank you very much!!! This solved my problem. This is very helpful! 
> Best regards,
> Yuesu
>
> On Wed, Aug 5, 2020 at 9:33 PM Bruno Blais  > wrote:
>
>> The best way to proceed is to set the boundary conditions using Physical 
>> Entities in GMSH.
>> An example is always better. For example, the following GMSH code :
>>
>>  lc = 2.0e-1;
>>  lf = 2.0e-1;
>>  RO=1;
>>  RI=0.25;
>>  
>>  Point(0) = {0, 0, 0, lc};
>>  Point(1) = {RO, 0, 0, lc};
>>  Point(2) = {0, -RO , 0, lc};
>>  Point(3) = {-RO, 0, 0, lc};
>>  Point(4) = {0, RO, 0, lc};
>>  
>>  Point(5) = {RI, 0, 0, lf};
>>  Point(6) = {0, -RI , 0, lf};
>>  Point(7) = {-RI, 0, 0, lf};
>>  Point(8) = {0, RI, 0, lf};
>>  
>>  Circle(1)={1,0,2};
>>  Circle(2)={2,0,3};
>>  Circle(3)={3,0,4};
>>  Circle(4)={4,0,1};
>>  
>>  Circle(5)={5,0,6};
>>  Circle(6)={6,0,7};
>>  Circle(7)={7,0,8};
>>  Circle(8)={8,0,5};
>>  
>>  Line Loop(1) = {1,2,3,4};
>>  Line Loop(2) = {5,6,7,8};
>>  
>>  Plane Surface(1) = {1,2} ;
>>  //Transfinite Surface{1}={1,2,3,4};
>>  Recombine Surface{1,2};
>>  
>>  // Creates a physical entity 1 (i.e. for a BC)
>>  //Physical Point(1) = {1,2} ;
>>  Physical Line(0)={1,2,3,4};
>>  Physical Line(1)={5,6,7,8};
>>  
>>  Physical Surface(0) = {1};
>>
>>
>> Produces a 2D mesh of a single surface with id 0. It has two boundary 
>> conditions (0) and (1). One on the outer circle and the other on the inner 
>> circle. We use it to simulate a Taylor-Couette flow.
>> If you wanted to use a 3D mesh, then you would make a Physical Volume(0) 
>> and you would use the Physical Surface for the boundary conditions
>>
>> Best
>> Bruno
>>
>>
>>
>> On Wednesday, 5 August 2020 20:36:09 UTC-4, yuesu jin wrote:
>>>
>>> Dear all,
>>>   I want to set up one Dirchlet boundary condition and one Neumann 
>>> boundary condition on a 2d  mesh, which is generated by gmsh. How can I set 
>>> up the mesh file in order to tell dealii with the boundary id when I input 
>>> the .msh file into dealii? Thanks!
>>> Best regards,
>>> -- 
>>> Yuesu Jin,
>>> Ph.D student,
>>> University of Houston,
>>> College of Natural Sciences and Mathematics,
>>> Department of Earth and Atmospheric Sciences,
>>> Houston, Texas 77204-5008 
>>> 346-404-2062
>>>
>>> -- 
>> 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/756546c6-95a9-4b71-a58b-772a0884bb9do%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/756546c6-95a9-4b71-a58b-772a0884bb9do%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
>
>
> -- 
> Yuesu Jin,
> Ph.D student,
> University of Houston,
> College of Natural Sciences and Mathematics,
> Department of Earth and Atmospheric Sciences,
> Houston, Texas 77204-5008 
> 346-404-2062
>
>

-- 
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/37e3ef06-2abe-49c7-b172-b23222025a94o%40googlegroups.com.

[deal.II] Re: Set up boundary id in gmsh

[Bruno Blais]

The best way to proceed is to set the boundary conditions using Physical 
Entities in GMSH.
An example is always better. For example, the following GMSH code :

 lc = 2.0e-1;
 lf = 2.0e-1;
 RO=1;
 RI=0.25;
 
 Point(0) = {0, 0, 0, lc};
 Point(1) = {RO, 0, 0, lc};
 Point(2) = {0, -RO , 0, lc};
 Point(3) = {-RO, 0, 0, lc};
 Point(4) = {0, RO, 0, lc};
 
 Point(5) = {RI, 0, 0, lf};
 Point(6) = {0, -RI , 0, lf};
 Point(7) = {-RI, 0, 0, lf};
 Point(8) = {0, RI, 0, lf};
 
 Circle(1)={1,0,2};
 Circle(2)={2,0,3};
 Circle(3)={3,0,4};
 Circle(4)={4,0,1};
 
 Circle(5)={5,0,6};
 Circle(6)={6,0,7};
 Circle(7)={7,0,8};
 Circle(8)={8,0,5};
 
 Line Loop(1) = {1,2,3,4};
 Line Loop(2) = {5,6,7,8};
 
 Plane Surface(1) = {1,2} ;
 //Transfinite Surface{1}={1,2,3,4};
 Recombine Surface{1,2};
 
 // Creates a physical entity 1 (i.e. for a BC)
 //Physical Point(1) = {1,2} ;
 Physical Line(0)={1,2,3,4};
 Physical Line(1)={5,6,7,8};
 
 Physical Surface(0) = {1};


Produces a 2D mesh of a single surface with id 0. It has two boundary 
conditions (0) and (1). One on the outer circle and the other on the inner 
circle. We use it to simulate a Taylor-Couette flow.
If you wanted to use a 3D mesh, then you would make a Physical Volume(0) 
and you would use the Physical Surface for the boundary conditions

Best
Bruno



On Wednesday, 5 August 2020 20:36:09 UTC-4, yuesu jin wrote:
>
> Dear all,
>   I want to set up one Dirchlet boundary condition and one Neumann 
> boundary condition on a 2d  mesh, which is generated by gmsh. How can I set 
> up the mesh file in order to tell dealii with the boundary id when I input 
> the .msh file into dealii? Thanks!
> Best regards,
> -- 
> Yuesu Jin,
> Ph.D student,
> University of Houston,
> College of Natural Sciences and Mathematics,
> Department of Earth and Atmospheric Sciences,
> Houston, Texas 77204-5008 
> 346-404-2062
>
>

-- 
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/756546c6-95a9-4b71-a58b-772a0884bb9do%40googlegroups.com.

[deal.II] Re: Getting started and posting guidelines for new users

[Bruno Blais]

Dear Kaleem,
It would be good to open your own topic with your question so that you may 
get help.
Best
Bruno


On Wednesday, 5 August 2020 05:30:44 UTC-4, kaleem iqbal wrote:
>
> The following error shown in step-49 during make run.
> Exception on processing: 
>
> 
> An error occurred in line <1430> of file 
>  
> in function
> void dealii::GridIn::read_msh(std::istream&) [with int 
> dim = 2; int spacedim = 2; std::istream = std::basic_istream]
> The violated condition was: 
> in
> Additional information: 
> An input/output error has occurred. There are a number of reasons why 
> this may be happening, both for reading and writing operations.
>
> If this happens during an operation that tries to read data: First, you 
> may be trying to read from a file that doesn't exist or that is not 
> readable given its file permissions. Second, deal.II uses this error at 
> times if it tries to read information from a file but where the information 
> in the file does not correspond to the expected format. An example would be 
> a truncated file, or a mesh file that contains not only sections that 
> describe the vertices and cells, but also sections for additional data that 
> deal.II does not understand.
>
> If this happens during an operation that tries to write data: you may be 
> trying to write to a file to which file or directory permissions do not 
> allow you to write. A typical example is where you specify an output file 
> in a directory that does not exist.
> 
>
> Aborting!
> 
> CMakeFiles/run.dir/build.make:57: recipe for target 'CMakeFiles/run' failed
> make[3]: *** [CMakeFiles/run] Error 1
> CMakeFiles/Makefile2:131: recipe for target 'CMakeFiles/run.dir/all' failed
> make[2]: *** [CMakeFiles/run.dir/all] Error 2
> CMakeFiles/Makefile2:138: recipe for target 'CMakeFiles/run.dir/rule' 
> failed
> make[1]: *** [CMakeFiles/run.dir/rule] Error 2
> Makefile:144: recipe for target 'run' failed
> make: *** [run] Error 2
>
> Regard's
> Kaleem iqbal
>

-- 
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/7cce5820-4ae3-469e-9a03-749805fae759o%40googlegroups.com.

[deal.II] Re: New docker image for deal.II 9.2.0

[Bruno Blais]

Awesome!
Thank you so much for taking the time to prepare these images Luca. They 
are amazingly helpful for those using Travis_CI for their deal.II derived 
codes :).


On Wednesday, 15 July 2020 05:04:31 UTC-4, luca.heltai wrote:
>
> Dear all, 
>
> I’d like to inform you that we now have also docker images with deal.II 
> 9.2.0, based on ubuntu bionic (18.04.1), and based on ubuntu focal 
> (20.04.1). 
>
> You can find them with 
>
> docker pull dealii/dealii:v9.2.0-bionic 
> docker pull dealii/dealii:v9.2.0-focal  (dealii/dealii:latest) 
>
> Best, 
> Luca. 
>
>

-- 
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/77ba1dd9-9623-4392-9398-f33631a0328do%40googlegroups.com.

Re: [deal.II] Interest in a step-12-like DG tutorial but for the advection-diffusion equation?

[Bruno Blais]

Dear Timo,
A step-12-like version of Step-39 for the Laplace equation would be very 
interesting. I would be glad to help on making one of those. The same could 
be said about the Stokes equations.

Ideally, once that is done, i would like to try and make one for stabilized 
Navier-Stokes for both CG and DG, but I feel I need a bit more of 
experience with writing steps before I move on to that.

So if you feel you need any help with the two aforementioned steps, I would 
be glad to help. I really like what you have done with the 
FEInterfaceValues, so if we can make it more popular, why not ? :)!
I am also integrating teaching DG within one of my graduate classes and I 
feel using deal.II for that will be really helpful.

Best
Bruno




On Friday, 10 July 2020 07:39:08 UTC-4, Timo Heister wrote:
>
> Hi Bruno, 
>
> I introduced FEInterfaceValues because I was frustrated by the lack of 
> DG examples and as a consequence DG users in the community and I 
> realized that a more intuitive way to teach it was necessary. Of 
> course we need more than 1 or 2 examples to do that. 
> I have a lot of unfinished stuff on my machine including working 
> Laplace and Stokes examples (see my todo list here 
> https://github.com/dealii/dealii/issues/8884 ) that I haven't gotten 
> back to yet. 
>
> Let me know if you want to work together on something or if you have 
> any ideas or suggestions. 
>
> Best, 
> Timo 
>
>
> On Thu, Jul 9, 2020 at 8:50 PM Bruno Blais  > wrote: 
> > 
> > Dear all, 
> > I hope you are well. 
> > I have been playing with the DG methods within deal.II lately (which has 
> been lots of fun). 
> > I really enjoy the way step-12 is written, especially because of the use 
> of the FEInterfaceValue class. In my opinion it makes low level stuff very 
> easily accessible, yet it is relatively easy to understand. 
> > However, I am under the impression that there are no similar test for 
> elliptic equations in which you need to use Nitsche method + Symmetric 
> Interior Penalty for the boundary + inner faces respectively. I know of 
> step-39, but it does not use the same "low-levelish" features. 
> > I have written such a code for my own pleasure for the laplace equation 
> and I am currently working on its advection-diffusion version (mostly 
> following Larson Chapter 14). 
> > 
> > Would it be interesting to turn such a case into a deal.II step (that 
> would follow in step-12 footstep) or does the community feel that the 
> existing steps for DG sufficiently cover the ground? 
> > 
> > Best 
> > Bruno 
> > 
> > -- 
> > 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/b931b23c-b284-4047-a887-d500ac9e8486o%40googlegroups.com.
>  
>
>
>
>
> -- 
> Timo Heister 
> http://www.math.clemson.edu/~heister/ 
>

-- 
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/9340f75c-a631-46e0-ac2e-400c6d14ff57o%40googlegroups.com.

[deal.II] Interest in a step-12-like DG tutorial but for the advection-diffusion equation?

[Bruno Blais]

Dear all,
I hope you are well.
I have been playing with the DG methods within deal.II lately (which has 
been lots of fun).
I really enjoy the way step-12 is written, especially because of the use of 
the FEInterfaceValue class. In my opinion it makes low level stuff very 
easily accessible, yet it is relatively easy to understand.
However, I am under the impression that there are no similar test for 
elliptic equations in which you need to use Nitsche method + Symmetric 
Interior Penalty for the boundary + inner faces respectively. I know of 
step-39, but it does not use the same "low-levelish" features.
I have written such a code for my own pleasure for the laplace equation and 
I am currently working on its advection-diffusion version (mostly following 
Larson Chapter 14).

Would it be interesting to turn such a case into a deal.II step (that would 
follow in step-12 footstep) or does the community feel that the existing 
steps for DG sufficiently cover the ground?

Best
Bruno

-- 
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/b931b23c-b284-4047-a887-d500ac9e8486o%40googlegroups.com.

[deal.II] Re: particle parallelization

[Bruno Blais]

Dear Shahab,
No, if you change a property of a particle on proc 0 and this particle is a 
ghost cell on proc 0, this property will not change on proc 1 (which would 
be its owner).
You would need to change the property on proc 1 and then I think the update 
on proc 0 will only be done when you call the exchange_ghost_particles 
function of the ParticleHandler class


On Sunday, 5 July 2020 23:58:16 UTC-4, Shahab Golshan wrote:
>
> Hello,
> I have a few questions on the parallelization aspects of a system 
> containing particles and I would appreciate it if someone could answer them.
> The first one is: Let's say I'm using 2 CPUs for a simulation. If I change 
> a property of a ghost particle (located in a ghost cell for CPU0), will the 
> property of the real particle on the other CPU change?
> Thanks in advance
> Shahab
>

-- 
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/ae17f8ab-560a-45b4-9fcc-0a1e5c0a0467o%40googlegroups.com.

[deal.II] Re: Xcode as an IDE for deal.II

[Bruno Blais]

Dear Siva,
I am unsure if I understand your question correctly.
When you run deal.II from within Xcode, the output of std::cout is printed 
to the console or it is not printed?


On Thursday, 11 June 2020 15:29:06 UTC-4, Siva Nadarajah wrote:
>
> Dear All, 
> My group has been using deal.ii for the past two years. I just switched to 
> using XCode. I have a problem where std::cout does write out to the console 
> within Xcode, but I do see it when the code runs on command line. Any 
> thoughts?
> Best Regards,
> Siva
>

-- 
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/80bc5391-9ea6-4f65-a07a-f87aed61fa36o%40googlegroups.com.

Re: [deal.II] Particle Class Implementation

[Bruno Blais]

Dear Victoria,

In additional to what Jean-Paul suggested, you can look at the preliminary 
version of step-68 which does exactly what you would like to achieve with 
particles.
The code is available on the following pull request. Rene and I have put 
some work into it and it works quite well in its current state:
https://github.com/dealii/dealii/pull/10308

Best
Bruno


On Tuesday, 9 June 2020 16:06:15 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Victoria,
>
> You’re on the right track, but it looks like you got the syntax a bit 
> wrong. To create a particle you need to write something like
>
> Particles::Particle particle;
>
> You’ll note that the specific constructor that is called is inferred from 
> the arguments given to the instance of the class that’s being created (in 
> this case, no arguments —> default constructor).
>
> There is a new tutorial, namely step-70 
> , that uses 
> particles. I don’t know whether or not it’s pitched at the right level for 
> you to get to grips with some of the fundamentals (at first glance, it 
> seems to involve a number of more advanced techniques). If it’s not quite 
> tractable for you, then you could also take a look at some of the particle 
> tests  in 
> the test suite, such as this one 
> ,
>  
> that might provide snippets of information that you need as you build up 
> your understanding of this particular feature of the library.
>
> I hope that helps you a bit.
> Best,
> Jean-Paul
>
> On 09 Jun 2020, at 17:14, Victoria W. > 
> wrote:
>
> Hello, 
>
> I'm a new deal.ii and C++ user, so this question might be a bit basic, but 
> I was wondering how my implementation of the built-in deal.ii particle 
> class should look.  Ultimately I'm looking to track a particle in laminar 
> flow.  So far I've tried creating a simple particle at the origin after 
> generating my geometry with the GridGenerator class.  Calling the particle 
> class in my code looks like this: 
>
> Particles::Particle::Particle() 
>
> as per the documentation on dealii.org.  However I get the error message: 
> error: dependent-name ‘dealii::Particles::Particle::Particle’ is 
> parsed as a non-type, but instantiation yields a type.  How should I be 
> implementing this class to construct a particle, what type declaration 
> should new particles be declared as, and are there other documentation or 
> resources I should be looking at?  
>
> 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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/43d0b22c-d736-464d-9d83-da4d819eca90o%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/579f09e4-b8ab-4197-8b05-dbd5e6fde291o%40googlegroups.com.

[deal.II] Strategy to snap the boundary of a triangulation to a manifold

[Bruno Blais]

Dear all,
I hope you are doing well.

In my endless quest for robust mesh generation of hex meshes using GMSH, I 
have managed to come up with a very robust strategy to generate hex-only 
meshes
My only issue (which is a major one) is that this implies that my 
decomposition from tet to hex adds nodes that are not "snapped" to the 
boundary, but that are only linear interpolation of the other node on the 
triangular faces.
Consequently, my quest remains unfulfilled.

Meshing through high-order and snapping the additional node to a high-order 
mesh from within GMSH is very troublesome and not very robust (and also 
very time consuming). However, an idea came to mind.
I was wondering if there could be an easy way to "snap" my faces to the 
manifold to which they belong.

My problem is thus the following:
- Given a triangulation and a manifold
- Some nodes are exactly on the manifolds (the original nodes of the tets) 
and some are not (the added nodes in the subdivision)
- What would be the best way to deform mesh so that the non-conforming node 
get deformed to the position which would be implied by the manifold? I 
think I could also make the process more robust by solving an additional 
elasticity equation during the deformation to deform the entire mesh 
instead of just the nodes close to the manifold.


Would any of you have a suggestion on how best to achieve the deformation 
of the nodes to match the manifold?


-- 
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/8d966870-4e33-4915-ad6c-2342018d82c3o%40googlegroups.com.

Re: [deal.II] Re: Almost installed Deal II, but there is a fly in the ointment!

[Bruno Blais]

Dear Wilmar,

This is not a problem. Members of the deal.II user group are very friendly 
and very helpful. I was helped tremendously by a lot of people (JP 
Pelteret, Wolfgang, Luca) when I started. 
Feel free to try and help others when you become more accustomed.

PS : If you must know, I am from the great cold northern country of Canada. 
Land of noble and funny animals such as the moose and the beaver.

On Sunday, 7 June 2020 19:52:33 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>
> Yes, I manage to install it!
>
> Thank you Bruno. 
>
> I don't know where you are from, but for sure it must be a very good place.
>
> Three months trying to install this very important software. It was 
> machines and operator limitations. It will be very important to my 
> Electrical Engineering graduation.
>
> Now I will resume The Finite Element Course on Cursera. 
>
> The step-1 example compiled and function!
>
> Sore for my naives questions. You are a nice person. Thank you again.
>
> Sincerely,
>
> Wilmar.
>
> Em dom., 7 de jun. de 2020 às 17:51, Wilmar Alves Cruvinel Lima <
> wilmar8...@gmail.com > escreveu:
>
>> Dear Bruno,
>>
>> I did only the commandmake. 
>>
>> It runs relatively fast. Now I put 
>>
>> make install
>>
>> After 1 1/2 hour it is on 70% of Building CXX object source...
>>
>> Until now noting is on /home/wilmar/dealii_library
>>
>> I hope that it will finish ok!
>>
>> Thank you again.
>>
>> Best regards.
>>
>> Wilmar. 
>>
>> Em Dom, 7 de jun de 2020 15:37, Bruno Blais > > escreveu:
>>
>>> Did you compile the library by using the "make install" command?
>>> Cmake only configures the build but does not compile the library per say.
>>> Best
>>> Bruno
>>>
>>>
>>> On Sunday, 7 June 2020 14:17:44 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>>>>
>>>> Dear Bruno,
>>>>
>>>> Great! It functioned...
>>>>
>>>> I put the commands:
>>>> wilmar@linuxmint:~/build$ *mkdir dealii_library*
>>>>
>>>> So, I put the command and it runned:
>>>> wilmar@linuxmint:~/build$ *cmake /home/wilmar/dealii-9.1.1 
>>>> -DCMAKE_INSTALL_PREFIX=/home/wilmar/dealii_library*
>>>> -- This is CMake 3.10.2
>>>> -- 
>>>> -- Include /home/wilmar/dealii-9.1.1/cmake/setup_external_macros.cmake
>>>> -- Include /home/wilmar/dealii-9.1.1/cmake/macros/macro_add_flags.cmake
>>>> -- Include 
>>>> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_compiler_setup.cmake
>>>> -- Include 
>>>> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_cxx_compiler_bug.cmake
>>>> -- Include 
>>>> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_mpi_interface.cmake
>>>> . . .
>>>>
>>>> At the end, nothing was on  */home/wilmar/dealii_library.  *Maybe it 
>>>> is on dealii-9.1.1 but I didn't localize. 
>>>>
>>>> In the last step I did something wrong?
>>>>
>>>> Em dom., 7 de jun. de 2020 às 13:41, Bruno Blais  
>>>> escreveu:
>>>>
>>>>> So there are two things you must keep in mind when you run cmake
>>>>>
>>>>>
>>>>> The command should be something similar to 
>>>>> cmake /path/to/dealii/sources 
>>>>> -DCMAKE_INSTALL_PREFIX=/path/to/where/you/want/to/install
>>>>>
>>>>> The first folder you point to is the folder you cloned the deal.II 
>>>>> sources to (which will contain a CMakeLists.txt
>>>>> The second folder you point to is where you want to install. Keep in 
>>>>> mind, the install directory should be different from the build directory.
>>>>> What happens is that you go from the sources to a build folder, and 
>>>>> then you install the library. The build folder is an intermediary folder.
>>>>>
>>>>> Best
>>>>> Bruno
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Sunday, 7 June 2020 12:06:38 UTC-4, Wilmar Alves Cruvinel Lima 
>>>>> wrote:
>>>>>>
>>>>>> Thank you Bruno,
>>>>>>
>>>>>> Yes, the valid path, I think, is  */home/wilmar/build . *Orientation 
>>>>>> on *https://www.dealii.org/current/readme.html 
>>>>>> <https://www.dealii.org/current/readme.html>* inform a different 
&

Re: [deal.II] Re: Almost installed Deal II, but there is a fly in the ointment!

[Bruno Blais]

Did you compile the library by using the "make install" command?
Cmake only configures the build but does not compile the library per say.
Best
Bruno


On Sunday, 7 June 2020 14:17:44 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>
> Dear Bruno,
>
> Great! It functioned...
>
> I put the commands:
> wilmar@linuxmint:~/build$ *mkdir dealii_library*
>
> So, I put the command and it runned:
> wilmar@linuxmint:~/build$ *cmake /home/wilmar/dealii-9.1.1 
> -DCMAKE_INSTALL_PREFIX=/home/wilmar/dealii_library*
> -- This is CMake 3.10.2
> -- 
> -- Include /home/wilmar/dealii-9.1.1/cmake/setup_external_macros.cmake
> -- Include /home/wilmar/dealii-9.1.1/cmake/macros/macro_add_flags.cmake
> -- Include 
> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_compiler_setup.cmake
> -- Include 
> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_cxx_compiler_bug.cmake
> -- Include 
> /home/wilmar/dealii-9.1.1/cmake/macros/macro_check_mpi_interface.cmake
> . . .
>
> At the end, nothing was on  */home/wilmar/dealii_library.  *Maybe it is 
> on dealii-9.1.1 but I didn't localize. 
>
> In the last step I did something wrong?
>
> Em dom., 7 de jun. de 2020 às 13:41, Bruno Blais  > escreveu:
>
>> So there are two things you must keep in mind when you run cmake
>>
>>
>> The command should be something similar to 
>> cmake /path/to/dealii/sources 
>> -DCMAKE_INSTALL_PREFIX=/path/to/where/you/want/to/install
>>
>> The first folder you point to is the folder you cloned the deal.II 
>> sources to (which will contain a CMakeLists.txt
>> The second folder you point to is where you want to install. Keep in 
>> mind, the install directory should be different from the build directory.
>> What happens is that you go from the sources to a build folder, and then 
>> you install the library. The build folder is an intermediary folder.
>>
>> Best
>> Bruno
>>
>>
>>
>>
>> On Sunday, 7 June 2020 12:06:38 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>>>
>>> Thank you Bruno,
>>>
>>> Yes, the valid path, I think, is  */home/wilmar/build . *Orientation on 
>>> *https://www.dealii.org/current/readme.html 
>>> <https://www.dealii.org/current/readme.html>* inform a different 
>>> directory from unpacked Deal.II.
>>>
>>> So, I returned to run cmake again, with the command:
>>> wilmar@linuxmint:~build$ 
>>> *cmake -DCMAKE_INSTALL_PREFIX=/home/wilmar/build*
>>> but cmake missed the CMakeLists.txt:   
>>> *CMake Error: The source directory "/home/wilmar/build" does not appear 
>>> to contain CMakeLists.txt*
>>>
>>> Yes, CMakeLists.txt isn't on /home/wilmar/build. I don't know where the 
>>> CMakeLists.txt for this compilation is. 
>>>
>>> Thank you in advance for orientation!
>>>
>>> Wilmar (Brazil)  -  Running from Covid!
>>>
>>>
>>> Em dom., 7 de jun. de 2020 às 12:03, Bruno Blais  
>>> escreveu:
>>>
>>>> Hello,
>>>> You are getting this error because the installation path that you 
>>>> specified is not valid.
>>>> When you did you cmake command, you specific the installation path to 
>>>> :  */path/to/install/dir/share/deal.II/scripts*
>>>> This path on your linux machine is not a valid path, and consequently, 
>>>> this folder cannot be created.
>>>> You need to run cmake again and specify a valid installation path 
>>>> (something like /home/your_user/dealii/ for instance).
>>>> You are very close to getting it to fully compile :)
>>>>
>>>> Best
>>>> Bruno
>>>>
>>>>
>>>> On Sunday, 7 June 2020 09:54:39 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>>>>>
>>>>> I tried to install Deal II (dealii) on a PC with 2 GB RAM, but after 
>>>>> the command:
>>>>>
>>>>> wilmar@linuxmint:~build$ *make --jobs=1 install  *(jobs=4 was bad)  
>>>>> it stops on [ 31%] Building CXX object 
>>>>> source/fe/CMakeFiles/obj_fe_release.dir/fe_q.cc.o.
>>>>>
>>>>> So, I tried on another PC with 4 GB RAM. It builds CXX object until 
>>>>> 100%:
>>>>>
>>>>> ...
>>>>> [100%] Building CXX object 
>>>>> source/non_matching/CMakeFiles/obj_non_matching_debug.dir/immersed_surface_quadrature.cc.o
>>>>> [100%] Built target obj_non_matching_debug
>>>>> Scanning dependencies of target deal_II.g
>>

Re: [deal.II] Re: Almost installed Deal II, but there is a fly in the ointment!

[Bruno Blais]

So there are two things you must keep in mind when you run cmake


The command should be something similar to 
cmake /path/to/dealii/sources 
-DCMAKE_INSTALL_PREFIX=/path/to/where/you/want/to/install

The first folder you point to is the folder you cloned the deal.II sources 
to (which will contain a CMakeLists.txt
The second folder you point to is where you want to install. Keep in mind, 
the install directory should be different from the build directory.
What happens is that you go from the sources to a build folder, and then 
you install the library. The build folder is an intermediary folder.

Best
Bruno




On Sunday, 7 June 2020 12:06:38 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>
> Thank you Bruno,
>
> Yes, the valid path, I think, is  */home/wilmar/build . *Orientation on 
> *https://www.dealii.org/current/readme.html 
> <https://www.dealii.org/current/readme.html>* inform a different 
> directory from unpacked Deal.II.
>
> So, I returned to run cmake again, with the command:
> wilmar@linuxmint:~build$ 
> *cmake -DCMAKE_INSTALL_PREFIX=/home/wilmar/build*
> but cmake missed the CMakeLists.txt:   
> *CMake Error: The source directory "/home/wilmar/build" does not appear to 
> contain CMakeLists.txt*
>
> Yes, CMakeLists.txt isn't on /home/wilmar/build. I don't know where the 
> CMakeLists.txt for this compilation is. 
>
> Thank you in advance for orientation!
>
> Wilmar (Brazil)  -  Running from Covid!
>
>
> Em dom., 7 de jun. de 2020 às 12:03, Bruno Blais  > escreveu:
>
>> Hello,
>> You are getting this error because the installation path that you 
>> specified is not valid.
>> When you did you cmake command, you specific the installation path to :  
>> */path/to/install/dir/share/deal.II/scripts*
>> This path on your linux machine is not a valid path, and consequently, 
>> this folder cannot be created.
>> You need to run cmake again and specify a valid installation path 
>> (something like /home/your_user/dealii/ for instance).
>> You are very close to getting it to fully compile :)
>>
>> Best
>> Bruno
>>
>>
>> On Sunday, 7 June 2020 09:54:39 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>>>
>>> I tried to install Deal II (dealii) on a PC with 2 GB RAM, but after the 
>>> command:
>>>
>>> wilmar@linuxmint:~build$ *make --jobs=1 install  *(jobs=4 was bad)  it 
>>> stops on [ 31%] Building CXX object 
>>> source/fe/CMakeFiles/obj_fe_release.dir/fe_q.cc.o.
>>>
>>> So, I tried on another PC with 4 GB RAM. It builds CXX object until 100%:
>>>
>>> ...
>>> [100%] Building CXX object 
>>> source/non_matching/CMakeFiles/obj_non_matching_debug.dir/immersed_surface_quadrature.cc.o
>>> [100%] Built target obj_non_matching_debug
>>> Scanning dependencies of target deal_II.g
>>> [100%] Building CXX object source/CMakeFiles/deal_II.g.dir/dummy.cc.o
>>> [100%] Linking CXX shared library ../lib/libdeal_II.g.so
>>> [100%] Built target deal_II.g
>>> Scanning dependencies of target obj_integrators_debug
>>> [100%] Built target obj_integrators_debug
>>> Scanning dependencies of target obj_integrators_release
>>> [100%] Built target obj_integrators_release
>>> Scanning dependencies of target obj_rol_inst
>>> [100%] Built target obj_rol_inst
>>> Scanning dependencies of target obj_rol_release
>>> [100%] Built target obj_rol_release
>>> Scanning dependencies of target obj_rol_debug
>>> [100%] Built target obj_rol_debug
>>>
>>> but, there is a fly in the ointment. I proceeds:
>>>
>>> Install the project...
>>> -- Install configuration: "DebugRelease"
>>> CMake Error at cmake/scripts/cmake_install.cmake:41 (file):
>>>   file cannot create directory: 
>>> /path/to/install/dir/share/deal.II/scripts.
>>>   Maybe need administrative privileges.
>>> Call Stack (most recent call first):
>>>   cmake_install.cmake:42 (include)
>>>
>>> Makefile:73: recipe for target 'install' failed
>>> make: *** [install] Error 1
>>>
>>> The command did not manage to create the directoty  
>>> */path/to/install/dir/share/deal.II/scripts*
>>>
>>> I feel that I crossed the ocean and died on the beach!   but the cause 
>>> wasn't n-Cov19!
>>>
>>> Any help?
>>>
>> -- 
>> 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 subs

[deal.II] Re: Almost installed Deal II, but there is a fly in the ointment!

[Bruno Blais]

Hello,
You are getting this error because the installation path that you specified 
is not valid.
When you did you cmake command, you specific the installation path to :  
*/path/to/install/dir/share/deal.II/scripts*
This path on your linux machine is not a valid path, and consequently, this 
folder cannot be created.
You need to run cmake again and specify a valid installation path 
(something like /home/your_user/dealii/ for instance).
You are very close to getting it to fully compile :)

Best
Bruno


On Sunday, 7 June 2020 09:54:39 UTC-4, Wilmar Alves Cruvinel Lima wrote:
>
> I tried to install Deal II (dealii) on a PC with 2 GB RAM, but after the 
> command:
>
> wilmar@linuxmint:~build$ *make --jobs=1 install  *(jobs=4 was bad)  it 
> stops on [ 31%] Building CXX object 
> source/fe/CMakeFiles/obj_fe_release.dir/fe_q.cc.o.
>
> So, I tried on another PC with 4 GB RAM. It builds CXX object until 100%:
>
> ...
> [100%] Building CXX object 
> source/non_matching/CMakeFiles/obj_non_matching_debug.dir/immersed_surface_quadrature.cc.o
> [100%] Built target obj_non_matching_debug
> Scanning dependencies of target deal_II.g
> [100%] Building CXX object source/CMakeFiles/deal_II.g.dir/dummy.cc.o
> [100%] Linking CXX shared library ../lib/libdeal_II.g.so
> [100%] Built target deal_II.g
> Scanning dependencies of target obj_integrators_debug
> [100%] Built target obj_integrators_debug
> Scanning dependencies of target obj_integrators_release
> [100%] Built target obj_integrators_release
> Scanning dependencies of target obj_rol_inst
> [100%] Built target obj_rol_inst
> Scanning dependencies of target obj_rol_release
> [100%] Built target obj_rol_release
> Scanning dependencies of target obj_rol_debug
> [100%] Built target obj_rol_debug
>
> but, there is a fly in the ointment. I proceeds:
>
> Install the project...
> -- Install configuration: "DebugRelease"
> CMake Error at cmake/scripts/cmake_install.cmake:41 (file):
>   file cannot create directory: /path/to/install/dir/share/deal.II/scripts.
>   Maybe need administrative privileges.
> Call Stack (most recent call first):
>   cmake_install.cmake:42 (include)
>
> Makefile:73: recipe for target 'install' failed
> make: *** [install] Error 1
>
> The command did not manage to create the directoty  
> */path/to/install/dir/share/deal.II/scripts*
>
> I feel that I crossed the ocean and died on the beach!   but the cause 
> wasn't n-Cov19!
>
> Any help?
>

-- 
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/cdc2ec84-0007-485b-9a8f-71daec1301ebo%40googlegroups.com.

[deal.II] Re: About metallic ductile damage

[Bruno Blais]

Dear Itachi,
Deal.II is significantly different than abaqus in the fact that you easily 
write you own equation and solve from it using whatever mean you find 
necessary.
However, it does not come readily built with a library of physical solvers 
or constitutive laws. Consequently, you will have to build your own weak 
form for your problem of interest (although this can be made much easier 
with automatic differentiation facilities)..

So consequently, any user material can be defined in deal.II, but it is 
significantly different than what you would obtained using a UMAT file for 
instance.

Best
Bruno


On Saturday, 6 June 2020 16:08:47 UTC-4, itachi Ezio wrote:
>
> Hi,
>
> I'm Abaqus user and fresh new about using deal.II. I'd like to ask a few 
> questions base on deal.II functions.
>
> I'm looking for the simulation for metallic ductile damage behavior with 
> deal.II ( for example with triaxiality and lode angle), is it possible?
>
> And how is it like to write user materials in deal.II, especially in 
> plasticity behaviors? Is it also quit sophisticated as in UMAT in Abaqus?
>
> Looking forward to the answers, thanks for that!  
>
>
>

-- 
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/c535b46f-c2a6-4db9-ac58-c25bca4dbbeco%40googlegroups.com.

Re: [deal.II] Re: Finding vertices on boundaries

[Bruno Blais]

I think David solution is more along the lines of Shahab original question 
(at least from my understanding) :).


On Friday, 1 May 2020 09:46:58 UTC-4, David Wells wrote:
>
> Hi Shahab,
>
> You are right that vertices themselves don't store any information - this 
> is an implementation detail of deal.II (points do not store any topological 
> information in a Triangulation, unlike lines, quadrilaterals, or hexahedra).
>
> Just to check: when you write 'vertices', are you looking for information 
> about the degrees of freedom on the boundary, or are you interested in the 
> getting a list of the Point objects that are on the boundary? 
> What Bruno describes will do the former. If you want to do the later then 
> something like this will work:
>
> Triangulation triangulation;
>
> // do something to set up the triangulation
>
> std::set> boundary_vertices;
> for (const auto  : triangulation.active_face_iterators())
>   if (face->at_boundary())
> for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v)
>   boundary_vertices.insert(face->vertex(v));
>
> You could also store Point *> objects if you want to manipulate 
> the points once you have all of them.
>
> Does this answer your question? Please post again if you still need help 
> figuring this out!
> 
> Best,
> David
>
>
> On Thu, Apr 30, 2020 at 11:28 PM Bruno Blais  > wrote:
>
>> Dear Shahab,
>> I think the best solution (I might be wrong) is to loop through the 
>> boundary faces and acquire the DOF indices that lie on a boundary face. You 
>> could then store them in an std::map or a similar structure that prevents 
>> data duplication.
>>
>>
>> You can create a vector and store the dof indices of a face in a 
>> following fashion that I believe should work.
>> Assuming you already have a vector such as:
>> std::vector face_dofs;
>>
>>
>>
>>
>> In the loops within the cell then the faces, assuming you are at a face 
>> *face* and that you have a FiniteElement fe, you can do:
>> face_dofs.resize(fe.dofs_per_face);
>> face->get_dof_indices(face_dofs, cell->active_fe_index());
>>
>> This will gather all the dofs indices of that given face. If that face is 
>> a boundary face it means you have gathered all of its dofs location.
>> I hope that points you in the right direction and helps you!
>> Best
>> Bruno
>>
>>
>>
>>
>> On Thursday, 30 April 2020 23:05:31 UTC-4, Shahab Golshan wrote:
>>>
>>> Dear all,
>>> I want to find all the vertices located on boundaries. It seems the 
>>> function at_boundary() does not work for vertices (points). Do you have any 
>>> suggestions how to obtain all the vertices located on the boundaries?
>>> Best,
>>> Shahab
>>>
>> -- 
>> 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/389defa8-9b6b-467e-a2f4-3d77abee31f0%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/389defa8-9b6b-467e-a2f4-3d77abee31f0%40googlegroups.com?utm_medium=email_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/ab8072ea-0d41-4577-a346-2bebfe4a7611%40googlegroups.com.

[deal.II] Re: Finding vertices on boundaries

[Bruno Blais]

Dear Shahab,
I think the best solution (I might be wrong) is to loop through the 
boundary faces and acquire the DOF indices that lie on a boundary face. You 
could then store them in an std::map or a similar structure that prevents 
data duplication.


You can create a vector and store the dof indices of a face in a following 
fashion that I believe should work.
Assuming you already have a vector such as:
std::vector face_dofs;




In the loops within the cell then the faces, assuming you are at a face 
*face* and that you have a FiniteElement fe, you can do:
face_dofs.resize(fe.dofs_per_face);
face->get_dof_indices(face_dofs, cell->active_fe_index());

This will gather all the dofs indices of that given face. If that face is a 
boundary face it means you have gathered all of its dofs location.
I hope that points you in the right direction and helps you!
Best
Bruno




On Thursday, 30 April 2020 23:05:31 UTC-4, Shahab Golshan wrote:
>
> Dear all,
> I want to find all the vertices located on boundaries. It seems the 
> function at_boundary() does not work for vertices (points). Do you have any 
> suggestions how to obtain all the vertices located on the boundaries?
> Best,
> Shahab
>

-- 
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/389defa8-9b6b-467e-a2f4-3d77abee31f0%40googlegroups.com.

[deal.II] Adding the manifold ID to a DataOutFaces

[Bruno Blais]

Dear all, 
I hope you are well.
I am currently working extensively with adaptive mesh refinement and I feel 
the need to output the 2D boundary surface of my 3D mesh (thus a dim=2, 
spacedim=3 type of problem)
I have realized there is an amazing Data out object that exists  
DataOutFaces which does exactly what I need.
I have figured out how to attach data to it and how to write parallel vtu 
files with it. However, I would like to add a data to this DataOutFace that 
would give out the manifold number (or the boundary id) of my boundary 
condition. However, since this DataOutFace only contains face of the 
surface, I am really unsure how to iterate on this structure.
I saw that there is a first_face() and next_face() function that should 
allow me to iterate over the face. However, this structure appears to be an 
std::pair and from this point I am confused as to how to use it.
Could somebody point me in the right direction? I have looked online, but 
it seems this class is not used extensively.

Thank you very much for your help,
Bruno

-- 
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/2effdbf6-a7e7-45f7-89cb-528865285460%40googlegroups.com.

Re: [deal.II] Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

Dear Jean-Paul, what you wrote in the FAQ is very clear. If I find 
additional elements as I play with the functionnality I will add it to the 
FAQ. Thanks for adding this section.
Best
Bruno


On Friday, 3 April 2020 16:11:05 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> I guess that it would be best to mention this in the tutorial itself, but 
> in the mean time I wrote a segment about it in the FAQ:
>
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-ensure-that-my-refined-grid-is-correct-when-using-cad-geometry
> Please feel free to edit it as you see fit. Never having used this 
> functionality in the library before, I hope that I have correctly 
> interpreted and massages what Luca had detailed.
>
> Best,
> Jean-Paul
>
> On 03 Apr 2020, at 15:55, Bruno Blais > 
> wrote:
>
> Dear Jean-Paul,
> Do you think it should be added to the wiki page or as additional 
> information to step-54? I have a bit on my hands, but i'd be glad to work 
> to make the information readily available. Clearly I would not have found 
> that out without Luca's help. 
> Best
> Bruno
>
>
> On Friday, 3 April 2020 09:43:33 UTC-4, Jean-Paul Pelteret wrote:
>>
>> Dear Bruno,
>>
>> I’m glad that you managed to solve this issue with Luca’s aid. I think 
>> that the information and workflow that Luca gave was really valuable. Would 
>> you care to add this to the Wiki page?
>>
>> Best,
>> Jean-Paul
>>
>> On 03 Apr 2020, at 06:03, Bruno Blais  wrote:
>>
>> Ciao Luca,
>> It works perfectly now even with regular GMSH triangulation (see movie). 
>> Now I understand the difference and how to set it correctly.
>> Thank you very much for the help :)!
>>
>>
>> On Thursday, 2 April 2020 11:31:42 UTC-4, luca.heltai wrote:
>>>
>>> Bruno, it seems like you are attaching your faces to the *boundary* 
>>> manifold. Let me try to explain what is happening. 
>>>
>>> In the code for step-54, there is a wire that is used to describe the 
>>> manifold of the *boundary* of the surface (a curve of dimension one 
>>> embedded in dimension three). This is generated by the wire that identifies 
>>> the boundary of the TopoDS shape of the face. The manifold attached to it 
>>> is an ArcLengthProjectionManifold, where mid points on the curve are added 
>>> by computing the distance in arc length, and taking the point in the middle 
>>> w.r.t. this length. 
>>>
>>> In your code with more than one cell, you are not assigning *just the 
>>> boundary faces* to the wire, but *all* faces (including the one in the 
>>> middle of your surface, delimitating the two cells). Intermediate points on 
>>> those lines will be projected (correctly) on the boundary wire (that’s what 
>>> your movie shows). 
>>>
>>> Of course any point that is in the middle of the surface does *not* 
>>> belong to the *boundary wire*, and you’ll get an exception when trying to 
>>> construct the midpoint of an *edge* that belongs to the interior of the 
>>> domain, with manifold id 2 (which corresponds to a manifold that describes 
>>> *the boundary* of your topology, not its surface). 
>>>
>>> With two cells, the internal edge (the only edge that does not belong to 
>>> the boundary) should *not* have id 2. If you set its id to 2, the result is 
>>> that a point which is in the middle of the edge, is actually in the midlle 
>>> of the two points *when running along the curve that describes the 
>>> boundary*, hence you get the distortion you show in the movie. 
>>>
>>> Try adding a check to the faces in which you set the manifold id to 2 
>>> only if the face is at the boundary. 
>>>
>>> The principle is: 
>>> 1. start from the lowest codimension objects, identify how to deform 
>>> them. In your case, cells of a Tria<2,3> are quads, that should deform as a 
>>> TopoDS_FACE. 
>>> 2. Attach your favorite manifold to the TopoDS_Face (I’d personally only 
>>> use NormalToMeshProjectionManifold) using 
>>> cell->set_all_manifold_ids(manifold_id) (Notice the use of 
>>> set_all_manifold_ids, and **not** set_manifold_id: you want all children of 
>>> these objects to belong to the same manifold, in particular you want all 
>>> faces that are internal, i.e., that are shared between two obects with the 
>>> same manifold id, to inherit the same manifold id). 
>>> 3. Go one codimension lower: in your case, curves (for 3d meshes, 
>>> surfaces). Follow the 

Re: [deal.II] Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

Dear Jean-Paul,
Do you think it should be added to the wiki page or as additional 
information to step-54? I have a bit on my hands, but i'd be glad to work 
to make the information readily available. Clearly I would not have found 
that out without Luca's help. 
Best
Bruno


On Friday, 3 April 2020 09:43:33 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> I’m glad that you managed to solve this issue with Luca’s aid. I think 
> that the information and workflow that Luca gave was really valuable. Would 
> you care to add this to the Wiki page?
>
> Best,
> Jean-Paul
>
> On 03 Apr 2020, at 06:03, Bruno Blais > 
> wrote:
>
> Ciao Luca,
> It works perfectly now even with regular GMSH triangulation (see movie). 
> Now I understand the difference and how to set it correctly.
> Thank you very much for the help :)!
>
>
> On Thursday, 2 April 2020 11:31:42 UTC-4, luca.heltai wrote:
>>
>> Bruno, it seems like you are attaching your faces to the *boundary* 
>> manifold. Let me try to explain what is happening. 
>>
>> In the code for step-54, there is a wire that is used to describe the 
>> manifold of the *boundary* of the surface (a curve of dimension one 
>> embedded in dimension three). This is generated by the wire that identifies 
>> the boundary of the TopoDS shape of the face. The manifold attached to it 
>> is an ArcLengthProjectionManifold, where mid points on the curve are added 
>> by computing the distance in arc length, and taking the point in the middle 
>> w.r.t. this length. 
>>
>> In your code with more than one cell, you are not assigning *just the 
>> boundary faces* to the wire, but *all* faces (including the one in the 
>> middle of your surface, delimitating the two cells). Intermediate points on 
>> those lines will be projected (correctly) on the boundary wire (that’s what 
>> your movie shows). 
>>
>> Of course any point that is in the middle of the surface does *not* 
>> belong to the *boundary wire*, and you’ll get an exception when trying to 
>> construct the midpoint of an *edge* that belongs to the interior of the 
>> domain, with manifold id 2 (which corresponds to a manifold that describes 
>> *the boundary* of your topology, not its surface). 
>>
>> With two cells, the internal edge (the only edge that does not belong to 
>> the boundary) should *not* have id 2. If you set its id to 2, the result is 
>> that a point which is in the middle of the edge, is actually in the midlle 
>> of the two points *when running along the curve that describes the 
>> boundary*, hence you get the distortion you show in the movie. 
>>
>> Try adding a check to the faces in which you set the manifold id to 2 
>> only if the face is at the boundary. 
>>
>> The principle is: 
>> 1. start from the lowest codimension objects, identify how to deform 
>> them. In your case, cells of a Tria<2,3> are quads, that should deform as a 
>> TopoDS_FACE. 
>> 2. Attach your favorite manifold to the TopoDS_Face (I’d personally only 
>> use NormalToMeshProjectionManifold) using 
>> cell->set_all_manifold_ids(manifold_id) (Notice the use of 
>> set_all_manifold_ids, and **not** set_manifold_id: you want all children of 
>> these objects to belong to the same manifold, in particular you want all 
>> faces that are internal, i.e., that are shared between two obects with the 
>> same manifold id, to inherit the same manifold id). 
>> 3. Go one codimension lower: in your case, curves (for 3d meshes, 
>> surfaces). Follow the same rule as above, setting all_manifold_ids on faces 
>> that you know should follow a known codimension one geometry (a known 
>> TopoDS_EDGE, or TopoDS_WIRE). In this case I’d only use 
>> ArcLengthProjectionManifold objects, attached to the wires that identify 
>> your geometry. 
>>
>> Doing things in this order guarantees that internal edges get the correct 
>> manifold id, which, in your case, is not happening. 
>>
>> Best, 
>> Luca. 
>>
>>
>>
>> > On 1 Apr 2020, at 17:32, Bruno Blais  wrote: 
>> > 
>> > So my investigation on this issue continues... 
>> > It seems the core of my issues stem from using a non-trivial mesh with 
>> more than one cell as a starting point. 
>> > if I start with a 2 cell spacedim=2 dim=3 mesh such as: 
>> > 
>> > 
>> >  
>> > 
>> > 
>> > 
>> > Which is identical to the second step of the adaptation. I get an 
>> aberrant result (see attached out_wrong.mp4). 
>> > 
>> > 
>> > 

[deal.II] Re: Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

Seems the image did not work,


On Wednesday, April 1, 2020 at 9:34:35 AM UTC-4, Bruno Blais wrote:
>
> So it seems part of my problem relate to the topology I was trying to 
> adapt OR to the staring complex mesh I was using.
> I have managed to make it work starting with a trivial mesh and using a 
> more simple IGES CAD. I will work on complexifying it from there. I will 
> try to keep the information here as it seems that feature has not been used 
> extensively :)!.
> What I have working now is this simple gif.
>
> [image: myimage.gif]
>
>
> On Tuesday, March 31, 2020 at 9:19:23 AM UTC-4, Bruno Blais wrote:
>>
>> Dear all,
>> I hope you are well in these uncertain times.
>>
>> I have been trying to use the OpenCASCADE library to set-up my manifolds. 
>> The goal of this is to be able to do CFD simulations in complex geometry, 
>> but use an initial very coarse mesh made with GMSH and use the dynamics 
>> mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. 
>> The geometry is complex, so fixing the manifold analytically appears to be 
>> a bad idea.
>> I started from step-54, which compiles and works well on my machines.
>> However, I am unable to "make it work" with my own generated Salome STEP 
>> and gmsh MSH files.
>> My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My 
>> mesh also contains numerous wires, so I disabled that aspect and only 
>> focused on refining with the faces which I label by looping over them.
>> However, no matter what occurs I get an error thrown:
>>
>>
>> 
>> An error occurred in line <114> of file 
>>  
>> in function
>> dealii::Point 
>> dealii::OpenCASCADE::NormalProjectionManifold> spacedim>::project_to_manifold(const dealii::ArrayView> dealii::Point >&, const dealii::Point&) const [with int 
>> dim = 2; int spacedim = 3]
>> The violated condition was: 
>> closest_point(sh, surrounding_points[i], tolerance) 
>> .distance(surrounding_points[i]) < std::max(tolerance * 
>> surrounding_points[i].norm(), tolerance)
>> Additional information: 
>> The point [ 0.1 0 0 ] is not on the manifold.
>>
>> Stacktrace:
>> ---
>> #0  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::OpenCASCADE::NormalProjectionManifold<2, 
>> 3>::project_to_manifold(dealii::ArrayView const, 
>> dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
>> #1  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::FlatManifold<2, 
>> 3>::get_new_point(dealii::ArrayView const, 
>> dealii::MemorySpace::Host> const&, dealii::ArrayView> dealii::MemorySpace::Host> const&) const
>> #2  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Manifold<2, 
>> 3>::get_new_point_on_line(dealii::TriaIterator> 3> > const&) const
>> #3  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> #4  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> #5  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
>> #6  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::DistortedCellList 
>> dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2,
>>  
>> 3>&, bool)
>> #7  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::execute_refinement()
>> #8  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
>> #9  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
>> dealii::Triangulation<2, 3>::refine_global(unsigned int)
>> #10  ./step-54: Step54::TriangulationOnCAD::refine_mesh()
>> #11  ./step-54: Step54::TriangulationOnCAD::run()
>> #12  ./step-54: main
>> 
>>
>> I have tried some elements :
>> - Playing with the tolerance of the CAD
>> - Writing my CAD from the OpenCASCADE STEP object and then opening it 
>> again in Salome to see that it is interpreted correctly +
>>
>> However, I always seem to face an issue where my CAD does not appear to 
>> generate a valid manifold. I have joined my 

[deal.II] Re: Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

So it seems part of my problem relate to the topology I was trying to adapt 
OR to the staring complex mesh I was using.
I have managed to make it work starting with a trivial mesh and using a 
more simple IGES CAD. I will work on complexifying it from there. I will 
try to keep the information here as it seems that feature has not been used 
extensively :)!.
What I have working now is this simple gif.

[image: myimage.gif]


On Tuesday, March 31, 2020 at 9:19:23 AM UTC-4, Bruno Blais wrote:
>
> Dear all,
> I hope you are well in these uncertain times.
>
> I have been trying to use the OpenCASCADE library to set-up my manifolds. 
> The goal of this is to be able to do CFD simulations in complex geometry, 
> but use an initial very coarse mesh made with GMSH and use the dynamics 
> mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. 
> The geometry is complex, so fixing the manifold analytically appears to be 
> a bad idea.
> I started from step-54, which compiles and works well on my machines.
> However, I am unable to "make it work" with my own generated Salome STEP 
> and gmsh MSH files.
> My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My 
> mesh also contains numerous wires, so I disabled that aspect and only 
> focused on refining with the faces which I label by looping over them.
> However, no matter what occurs I get an error thrown:
>
>
> 
> An error occurred in line <114> of file 
>  
> in function
> dealii::Point 
> dealii::OpenCASCADE::NormalProjectionManifold spacedim>::project_to_manifold(const dealii::ArrayView dealii::Point >&, const dealii::Point&) const [with int 
> dim = 2; int spacedim = 3]
> The violated condition was: 
> closest_point(sh, surrounding_points[i], tolerance) 
> .distance(surrounding_points[i]) < std::max(tolerance * 
> surrounding_points[i].norm(), tolerance)
> Additional information: 
> The point [ 0.1 0 0 ] is not on the manifold.
>
> Stacktrace:
> ---
> #0  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::OpenCASCADE::NormalProjectionManifold<2, 
> 3>::project_to_manifold(dealii::ArrayView const, 
> dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
> #1  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::FlatManifold<2, 
> 3>::get_new_point(dealii::ArrayView const, 
> dealii::MemorySpace::Host> const&, dealii::ArrayView dealii::MemorySpace::Host> const&) const
> #2  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Manifold<2, 
> 3>::get_new_point_on_line(dealii::TriaIterator 3> > const&) const
> #3  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> #4  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> #5  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
> #6  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::DistortedCellList 
> dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2,
>  
> 3>&, bool)
> #7  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::execute_refinement()
> #8  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
> #9  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::Triangulation<2, 3>::refine_global(unsigned int)
> #10  ./step-54: Step54::TriangulationOnCAD::refine_mesh()
> #11  ./step-54: Step54::TriangulationOnCAD::run()
> #12  ./step-54: main
> 
>
> I have tried some elements :
> - Playing with the tolerance of the CAD
> - Writing my CAD from the OpenCASCADE STEP object and then opening it 
> again in Salome to see that it is interpreted correctly +
>
> However, I always seem to face an issue where my CAD does not appear to 
> generate a valid manifold. I have joined my code as well as the salome HDF 
> file and gmsh GEO file used to generate the CAD and meshes respectively.
> I must admit I am quite at lost here as to why both don't coincide. Any 
> advices on where to proceed from there?
> Best :)!
> Bruno
>
>

-- 
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/a239784b-954a-4860-ab3a-acac05e4c709%40googlegroups.com.

[deal.II] Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54

[Bruno Blais]

Dear all,
I hope you are well in these uncertain times.

I have been trying to use the OpenCASCADE library to set-up my manifolds. 
The goal of this is to be able to do CFD simulations in complex geometry, 
but use an initial very coarse mesh made with GMSH and use the dynamics 
mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. 
The geometry is complex, so fixing the manifold analytically appears to be 
a bad idea.
I started from step-54, which compiles and works well on my machines.
However, I am unable to "make it work" with my own generated Salome STEP 
and gmsh MSH files.
My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My 
mesh also contains numerous wires, so I disabled that aspect and only 
focused on refining with the faces which I label by looping over them.
However, no matter what occurs I get an error thrown:



An error occurred in line <114> of file 
 
in function
dealii::Point 
dealii::OpenCASCADE::NormalProjectionManifold::project_to_manifold(const dealii::ArrayView >&, const dealii::Point&) const [with int 
dim = 2; int spacedim = 3]
The violated condition was: 
closest_point(sh, surrounding_points[i], tolerance) 
.distance(surrounding_points[i]) < std::max(tolerance * 
surrounding_points[i].norm(), tolerance)
Additional information: 
The point [ 0.1 0 0 ] is not on the manifold.

Stacktrace:
---
#0  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::OpenCASCADE::NormalProjectionManifold<2, 
3>::project_to_manifold(dealii::ArrayView const, 
dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
#1  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::FlatManifold<2, 
3>::get_new_point(dealii::ArrayView const, 
dealii::MemorySpace::Host> const&, dealii::ArrayView const&) const
#2  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Manifold<2, 
3>::get_new_point_on_line(dealii::TriaIterator > const&) const
#3  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
#4  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
#5  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
#6  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Triangulation<2, 3>::DistortedCellList 
dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2,
 
3>&, bool)
#7  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Triangulation<2, 3>::execute_refinement()
#8  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
#9  /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: 
dealii::Triangulation<2, 3>::refine_global(unsigned int)
#10  ./step-54: Step54::TriangulationOnCAD::refine_mesh()
#11  ./step-54: Step54::TriangulationOnCAD::run()
#12  ./step-54: main


I have tried some elements :
- Playing with the tolerance of the CAD
- Writing my CAD from the OpenCASCADE STEP object and then opening it again 
in Salome to see that it is interpreted correctly +

However, I always seem to face an issue where my CAD does not appear to 
generate a valid manifold. I have joined my code as well as the salome HDF 
file and gmsh GEO file used to generate the CAD and meshes respectively.
I must admit I am quite at lost here as to why both don't coincide. Any 
advices on where to proceed from there?
Best :)!
Bruno

-- 
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/fdc5899b-53d1-41ca-8efe-72742a52aed6%40googlegroups.com.


topology_adaptation.tar.gz
Description: Binary data

Re: [deal.II] How to visualise the simulation output in vtk format on a one dimensional mesh?

[Bruno Blais]

I concur with Wolfgang.
The same can be done with Paraview also. I have found that using the plot 
over line tool of paraview in 1D allows to obtain animation of 1D results 
in a very beautiful and user friendly fashion.
There is a Python VTK library by the way which allows easy (relatively) 
opening of vtk file and conversion into numpy arrays. It can be pretty 
handy, although the documentation is poor.

:)!


On Monday, 16 March 2020 13:00:57 UTC-4, Wolfgang Bangerth wrote:
>
> On 3/16/20 9:37 AM, Krishnakumar Gopalakrishnan wrote: 
> > 
> > Are we all agreeing that sophisticated VisIT and Paraview cannot produce 
> > simple 1D visualizations/animations? 
>
> No, we definitely do not agree :-) Attached is a visualization of the 
> output 
> when replacing all occurrences of <2> by <1> in step-3, using visit. I'm 
> showing the solution in color, elevated by the value of the solution. 
>
>
> > If so, there must be anAssert() 
> > pre-condition check for the DataOut::write_vtk() and 
> DataOut::write_vtu() 
> > functions of the library which return an exception with "dimension not 
> > implemented" error. 
>
> As others have pointed out, there are other reasons to output in these 
> file 
> formats that may not include visualization with Paraview only. 
>
> 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/cd00bc54-d8b9-43b6-928a-51e49e54ada9%40googlegroups.com.

Re: [deal.II] Re: solving stabilized Stokes

[Bruno Blais]

Dear Wolfgang,
I spent some time re-reading the theory and you are right, nothing shows 
that the convergence rate should be conserved when we have stabilized equal 
order elements. However it is interesting to note that when we use 
stabilization and we revert to LBB stable elements (Q2-Q1 for instance) the 
right order is recovered.

It seems that Richard was expecting that the order would be conserved 
however when PSPG stabilization is used. In our case, at least,I do not 
manage to reproduce that (and that is across the usage of two codes using 
GLS implementations).




On Wednesday, 11 March 2020 12:43:18 UTC-4, Wolfgang Bangerth wrote:
>
> On 3/11/20 8:54 AM, Bruno Blais wrote: 
> > Q1-Q1 using the GLS (SUPG +PSPG) stabilized solver of Lethe using a 
> > monolithic iterative solver and Q2-Q1 using the Grad-Div solver with a 
> > Schur completement solution strategy (similar to step-57 but using 
> > Trilinos). 
> > 
> > *Here are the results I obtain:* 
> > *Q1-Q1 - Velocity error and convergence* 
> >   cells   error 
> >   64 1.3282e-01- 
> >  256 3.4363e-02 1.95 
> > 1024 8.7362e-03 1.98 
> > 4096 2.1969e-03 1.99 
> >16384 5.5029e-04 2.00 
> >65536 1.3767e-04 2.00 
> >   262144 3.4426e-05 2.00 
> > 1048576 8.6075e-06 2.00 
> > /Conclusion : Order = p+1 recovered for the velocity!/ 
> > 
> > *Q1-Q1 - Pressure error and convergence* 
> > Refinement levelL2Error-PressureRatio 
> > 1   0.171975 
> > 2   0.0964683   1.33 
> > 3   0.0301739   1.78 
> > 4   0.00844618  1.89 
> > 5   0.0023758   1.885 
> > 6   0.000698421 1.84 
> > 7   0.000216927 1.79 
> > 8   7.07505e-05 1.75 
> > /Conclusion : Order = ??? The error does not reach the right asymptotic 
> > convergence rate./ 
>
> If I understand the chain of emails correctly, this is the question you 
> have: Using a Q1/Q1 element with GLS/PSPG, what convergence rates does 
> one expect. 
>
> You seem to suggest that it should be 2 for both velocity and pressure. 
> But is that true? What does the theory say? I would actually be quite 
> surprised if that's what you get for these kinds of problems. My gut 
> feeling is that you should expect 2=p+1 for the velocity, but something 
> between p and p+1 for the pressure. That would be consistent with other 
> stabilized approaches. 
>
> 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/9ef3def6-ce41-4787-9294-5e47c913b927%40googlegroups.com.

[deal.II] Re: solving stabilized Stokes

[Bruno Blais]

In addition to the above results. I ran the same MMS with the GLS 
stabilized solver but using Q2-Q1 and Q2-Q2.

What i obtain for Q2-Q1 is exactly what you would expect:
*Velocity*
Cell  error -
   64 1.4972e-02- 
  256 1.9438e-03 2.95 
 1024 2.4542e-04 2.99 
 4096 3.0755e-05 3.00 

*Pressure*
cell error -
64 0.0232558  -
256 0.00355176 2.55
1024 0.000831086 2.06
4096 0.000206113 2.01

And for Q2-Q2 the results become similar to Q1-Q1 in that the pressure 
convergence is lost:
Velocity
cells  error  
   64 1.5715e-02- 
  256 1.9762e-03 2.99 
 1024 2.4655e-04 3.00 
 4096 3.0792e-05 3.00 
16384 3.8480e-06 3.00

Pressure
cells  error  
   64 0.0796909 2.03
  256 0.0192026 2.01  
 1024 0.00476045 2.00
 4096 0.00118731 2.00
16384 0.000296642 2.00


On Wednesday, 11 March 2020 10:54:21 UTC-4, Bruno Blais wrote:
>
> Dear Richard,
> Thanks for your message! It is very interesting. Your results made me 
> doubt our own results, so I re-ran Error convergence analysis on a 
> manufactured solution ( an infinitely continuous one) where the domain is 
> bounded by purely Dirichlet boundary condition.
> I did the simulations with two approaches:
> Q1-Q1 using the GLS (SUPG +PSPG) stabilized solver of Lethe using a 
> monolithic iterative solver and Q2-Q1 using the Grad-Div solver with a 
> Schur completement solution strategy (similar to step-57 but using 
> Trilinos).
>
> *Here are the results I obtain:*
> *Q1-Q1 - Velocity error and convergence*
>  cells   error  
>  64 1.3282e-01- 
> 256 3.4363e-02 1.95 
>1024 8.7362e-03 1.98 
>4096 2.1969e-03 1.99 
>   16384 5.5029e-04 2.00 
>   65536 1.3767e-04 2.00 
>  262144 3.4426e-05 2.00 
> 1048576 8.6075e-06 2.00
> *Conclusion : Order = p+1 recovered for the velocity!*
>
> *Q1-Q1 - Pressure error and convergence*
> Refinement levelL2Error-PressureRatio
> 1   0.171975
> 2   0.0964683   1.33
> 3   0.0301739   1.78
> 4   0.00844618  1.89
> 5   0.0023758   1.885   
> 6   0.000698421 1.84
> 7   0.000216927 1.79
> 8   7.07505e-05 1.75
>
> *Conclusion : Order = ??? The error does not reach the right asymptotic 
> convergence rate. This is with a very low non-linear solver residual. 
> Clearly, There is a significant loss of accuracy in the pressure. The error 
> on the pressure is also orders of magnitudes higher than the grad-div 
> solver...*
>
>
> *Q2-Q1 - Velocity error and convergence*
> cells  error  
>64 1.4729e-02- 
>   256 1.9368e-03 2.93 
>  1024 2.4521e-04 2.98 
>  4096 3.0749e-05 3.00 
> 16384 3.8466e-06 3.00 
> 65536 4.8237e-07 3.00
> *Conclusion : Order = p+1 recovered for the velocity!*
>
> *Q2-Q1 - pressure error and convergence*
> Refinement levelL2Error-PressureRatio
> 1   0.0306665
> 2   0.003991442.77183454825005
> 3   0.000840952.178658165
> 4   0.0002063012.01899117611792
> 5   5.1481e-052.00182993535217
> 6   1.28777e-051.99942139684848
> *Conclusion : Order = p+1 recovered for the velocity!*
>
> If you want, I would be more than glad to discuss this issue even more. We 
> could organize maybe a skype call to discuss this? Clearly we are having 
> very very similar issues :)!
> Looking forward to more of this discussion.
> Best
> Bruno
>
>
>
> On Tuesday, 10 March 2020 12:05:26 UTC-4, Richard Schussnig wrote:
>>
>> Hi everyone,
>> I am back with good & bad news!
>> - Good news first, I implemented a parallel direct solver (mumps via 
>> petsc) and going from a 2x2 blocked system to 
>> a regular one (actually 1x1 still blocked system)
>> one just needs to adapt the setup phase, not reordering per component & 
>> not use the "get_view" 
>> but rather use all the same vector and matrix (types) with 1 block 
>> instead of 2. Once you realize that, it takes ~15 mins for a 
>> non-experienced deal.ii user like me.
>>
>> - Bad news is, that I still do not obtain the optimal convergence rate in 
>> the pressure, i.e., 2 for Q1Q1 PSPG-stabilized Stokes.
>>
>> The obtained rates in the Poisuille benchmark (quadratic velocity profile 
>> & linear pressure) using the direct solver are:
>>
>>  L2-errors in v & p: Q1Q1 + PSPG
>>
>> ---

[deal.II] Re: solving stabilized Stokes

[Bruno Blais]

Dear Richard,
Thanks for your message! It is very interesting. Your results made me doubt 
our own results, so I re-ran Error convergence analysis on a manufactured 
solution ( an infinitely continuous one) where the domain is bounded by 
purely Dirichlet boundary condition.
I did the simulations with two approaches:
Q1-Q1 using the GLS (SUPG +PSPG) stabilized solver of Lethe using a 
monolithic iterative solver and Q2-Q1 using the Grad-Div solver with a 
Schur completement solution strategy (similar to step-57 but using 
Trilinos).

*Here are the results I obtain:*
*Q1-Q1 - Velocity error and convergence*
 cells   error  
 64 1.3282e-01- 
256 3.4363e-02 1.95 
   1024 8.7362e-03 1.98 
   4096 2.1969e-03 1.99 
  16384 5.5029e-04 2.00 
  65536 1.3767e-04 2.00 
 262144 3.4426e-05 2.00 
1048576 8.6075e-06 2.00
*Conclusion : Order = p+1 recovered for the velocity!*

*Q1-Q1 - Pressure error and convergence*
Refinement levelL2Error-PressureRatio
1   0.171975
2   0.0964683   1.33
3   0.0301739   1.78
4   0.00844618  1.89
5   0.0023758   1.885   
6   0.000698421 1.84
7   0.000216927 1.79
8   7.07505e-05 1.75

*Conclusion : Order = ??? The error does not reach the right asymptotic 
convergence rate. This is with a very low non-linear solver residual. 
Clearly, There is a significant loss of accuracy in the pressure. The error 
on the pressure is also orders of magnitudes higher than the grad-div 
solver...*


*Q2-Q1 - Velocity error and convergence*
cells  error  
   64 1.4729e-02- 
  256 1.9368e-03 2.93 
 1024 2.4521e-04 2.98 
 4096 3.0749e-05 3.00 
16384 3.8466e-06 3.00 
65536 4.8237e-07 3.00
*Conclusion : Order = p+1 recovered for the velocity!*

*Q2-Q1 - pressure error and convergence*
Refinement levelL2Error-PressureRatio
1   0.0306665
2   0.003991442.77183454825005
3   0.000840952.178658165
4   0.0002063012.01899117611792
5   5.1481e-052.00182993535217
6   1.28777e-051.99942139684848
*Conclusion : Order = p+1 recovered for the velocity!*

If you want, I would be more than glad to discuss this issue even more. We 
could organize maybe a skype call to discuss this? Clearly we are having 
very very similar issues :)!
Looking forward to more of this discussion.
Best
Bruno



On Tuesday, 10 March 2020 12:05:26 UTC-4, Richard Schussnig wrote:
>
> Hi everyone,
> I am back with good & bad news!
> - Good news first, I implemented a parallel direct solver (mumps via 
> petsc) and going from a 2x2 blocked system to 
> a regular one (actually 1x1 still blocked system)
> one just needs to adapt the setup phase, not reordering per component & 
> not use the "get_view" 
> but rather use all the same vector and matrix (types) with 1 block instead 
> of 2. Once you realize that, it takes ~15 mins for a 
> non-experienced deal.ii user like me.
>
> - Bad news is, that I still do not obtain the optimal convergence rate in 
> the pressure, i.e., 2 for Q1Q1 PSPG-stabilized Stokes.
>
> The obtained rates in the Poisuille benchmark (quadratic velocity profile 
> & linear pressure) using the direct solver are:
>
>  L2-errors in v & p: Q1Q1 + PSPG
>
> --
> lvl   0: | e_v =1 | eoc_v =0 | e_p =1 
> | eoc_p =0 | 
>
> --
> lvl   1: | e_v = 0.531463 | eoc_v =  0.91196 | e_p = 0.479585 
> | eoc_p =  1.06014 | 
>
> --
> lvl   2: | e_v = 0.218196 | eoc_v =  1.28434 | e_p = 0.208925 
> | eoc_p =   1.1988 | 
>
> --
> lvl   3: | e_v =0.0661014 | eoc_v =  1.72287 | e_p =0.0674415 
> | eoc_p =  1.63128 | 
>
> --
> lvl   4: | e_v =0.0176683 | eoc_v =  1.90352 | e_p =0.0194532 
> | eoc_p =  1.79363 | 
>
> --
> lvl   5: | e_v =   0.00452588 | eoc_v =  1.96489 | e_p =   0.00549756 
> | eoc_p =  1.82315 | 
>
> --
> lvl   6: | e_v =   0.00114238 | eoc_v =  1.98616 | e_p =   0.00158448 
> | eoc_p =  1.79478 | 
>
> 

[deal.II] Re: solving stabilized Stokes

[Bruno Blais]

Dear Richard,
You are absolutely right in saying that if you have the Stokes equation 
your problem should be linear, even with the PSPG stabilization term.
Laplacian of u is only going to be zero for Q1-Q1 elements if your mesh is 
perfectly aligned with the axis of your domain. 
 However, this is what you mentioned so I do not think your error stems 
from there.
If it I were you, I would do as Wolfgang suggested and try to look at what 
happens when you use a direct solver to solve the system (say UMFPACK).
Which procedure are you using the define the "pressure constant" if you 
domain is enclosed?

Best
Bruno


On Wednesday, 19 February 2020 03:11:20 UTC-5, Richard Schussnig wrote:
>
> Hi Bruno,
> Thank you very much for taking the time to read my post & giving me some 
> ideas!
> I am already familiar with lethe, thats where i compared my implementation 
> of the residual to, but for Stokes, 
> we actually do not introduce a nonlinearity, since we do not have the 
> (nonlinear) convective term as in Navier Stokes, so it should actuallly be 
> fine, correct? 
> As far as I know, the Laplacian(u) in the residual does not lead to 
> nonlinearities, since the residual in PSPG for Stokes(!) is only r := 
> -nu*laplace(u) + grad(p) - f
> and we simply add on element level + int( r * grad(q) )dx
> Nonetheless, I solve it using Newton OR Picard+Aitken, giving me identical 
> results.
>
> That it takes longer to reach the asymptotic convergence I can confirm 
> (using my matlab code), 
> but using deal.ii and an iterative solver it just does not reach the point 
> (I tested up to 3 mio cells in 2d on just a rectangular grid with 
> poisuille flow ... didnt help)
> thanks also for that hint!
>
> Kind regards,
> Richard
>
> Am Mittwoch, 19. Februar 2020 02:49:26 UTC+1 schrieb Bruno Blais:
>>
>> Dear Richard,
>> We implement a quasi identical scheme in lethe (PSPG for continuity and 
>> SUPG for momentum). You can find the code on the github, it might help you :
>>
>> https://github.com/lethe-cfd/lethe/blob/master/source/solvers/gls_navier_stokes.cc
>>
>> I can confirm that you should be getting order p+1 for velocity. To be 
>> honest I have never really taken the time to look at pressure convergence, 
>> but I recall that it can take very fine mesh to reach the asymptotic 
>> convergence.
>>
>> Please keep in mind that the introduction of the PSPG term introduces a 
>> non-linearity in the solver except in the case where you rmesh is made of 
>> perfect rectangles that are aligned with the axis of the system.
>> Are you solving it using Newton's method?
>>
>> Best
>> Bruno
>>
>>
>> On Friday, 14 February 2020 07:34:09 UTC-5, Richard Schussnig wrote:
>>>
>>> Hi everyone!
>>> I am currently implementing some stationary Stokes solvers based on 
>>> step-55.
>>> Therein, Taylor-Hood elements are being used. One can check the optimal 
>>> order of 
>>> convergence easily comparing to the Kovasznay or Poisuille flow 
>>> solution, the first one being already implemented in this step.
>>> I went ahead and added PSPG and Bochev-Dohrmann stabilizations to check 
>>> the errors using those.
>>>
>>> For the latter, i get the expected suboptimal (!) convergence rates 
>>> (vel:2 and p~1.5-2).
>>> Using PSPG, one adds 
>>> tau * residual momentum dot gradient(q)
>>> (q being the pressure test function)
>>> per element like:
>>>
>>> local_matrix (i,j) += fe_values.JxW(q) * tau_sups * phi_grads_p [j] * 
>>> phi_grads_p [i];
>>>
>>> if (ADD_THE_LAPLACIAN)
>>> {
>>> local_matrix (i,j) += fe_values.JxW(q) * tau * (
>>>   - viscosity * phi_laplacians_v [j]
>>>   ) * phi_grads_p [i];
>>> }
>>>
>>> Using the laplacian of the velocity u defined as below:
>>>
>>> if (get_laplace_residual)
>>> {
>>> phi_hessians_v[k] = fe_hessians[velocities].hessian(k,q);
>>> for (int d = 0; d < dim; ++d)
>>> { phi_laplacians_v[k][d] = trace(phi_hessians_v[k][d]);
>>> }
>>>
>>> // Check laplacian, both give zero for rectangular grid and identical 
>>> values for distorted grids.
>>> unsigned int comp_k = fe.system_to_component_index(k).first;
>>>
>>> Tensor<1,dim> vector_laplace_v;
>>> vector_laplace_v = 0;
>>> if (comp_k>> {
>>> Tensor<2,dim> nonzero_hessian_k = 
>>> fe_hessians.shape_hessian_component(k,q,comp_k);
>>> vector_laplace_v [comp_k] += trace(nonzero_hessian_k);
>>> }
>

[deal.II] Re: solving stabilized Stokes

[Bruno Blais]

Dear Richard,
We implement a quasi identical scheme in lethe (PSPG for continuity and 
SUPG for momentum). You can find the code on the github, it might help you :
https://github.com/lethe-cfd/lethe/blob/master/source/solvers/gls_navier_stokes.cc

I can confirm that you should be getting order p+1 for velocity. To be 
honest I have never really taken the time to look at pressure convergence, 
but I recall that it can take very fine mesh to reach the asymptotic 
convergence.

Please keep in mind that the introduction of the PSPG term introduces a 
non-linearity in the solver except in the case where you rmesh is made of 
perfect rectangles that are aligned with the axis of the system.
Are you solving it using Newton's method?

Best
Bruno


On Friday, 14 February 2020 07:34:09 UTC-5, Richard Schussnig wrote:
>
> Hi everyone!
> I am currently implementing some stationary Stokes solvers based on 
> step-55.
> Therein, Taylor-Hood elements are being used. One can check the optimal 
> order of 
> convergence easily comparing to the Kovasznay or Poisuille flow solution, 
> the first one being already implemented in this step.
> I went ahead and added PSPG and Bochev-Dohrmann stabilizations to check 
> the errors using those.
>
> For the latter, i get the expected suboptimal (!) convergence rates (vel:2 
> and p~1.5-2).
> Using PSPG, one adds 
> tau * residual momentum dot gradient(q)
> (q being the pressure test function)
> per element like:
>
> local_matrix (i,j) += fe_values.JxW(q) * tau_sups * phi_grads_p [j] * 
> phi_grads_p [i];
>
> if (ADD_THE_LAPLACIAN)
> {
> local_matrix (i,j) += fe_values.JxW(q) * tau * (
>   - viscosity * phi_laplacians_v [j]
>   ) * phi_grads_p [i];
> }
>
> Using the laplacian of the velocity u defined as below:
>
> if (get_laplace_residual)
> {
> phi_hessians_v[k] = fe_hessians[velocities].hessian(k,q);
> for (int d = 0; d < dim; ++d)
> { phi_laplacians_v[k][d] = trace(phi_hessians_v[k][d]);
> }
>
> // Check laplacian, both give zero for rectangular grid and identical 
> values for distorted grids.
> unsigned int comp_k = fe.system_to_component_index(k).first;
>
> Tensor<1,dim> vector_laplace_v;
> vector_laplace_v = 0;
> if (comp_k {
> Tensor<2,dim> nonzero_hessian_k = 
> fe_hessians.shape_hessian_component(k,q,comp_k);
> vector_laplace_v [comp_k] += trace(nonzero_hessian_k);
> }
> Tensor<1,dim> diff;
> diff = phi_laplacians_v[k] - vector_laplace_v;
>
> if (diff.norm() > 1e-6 || vector_laplace_v.norm() > 1e-16)
> std::cout << "|divgrad_v| = " << std::setw(15) << vector_laplace_v.norm() 
> << " , diff = " << diff.norm() << std::endl;
>
> }
>
> Which also includes a second option to compute the wanted laplacian.
>
> Using a perfectly rectangular grid, the laplacian is actually zero for all 
> elements, thus everything reduces to just adding a pressure stiffness in 
> the pressure-pressure block.
> Nonetheless, i observe convergence rates of vel:2 and p:1.5-1.8 which is 
> suboptimal and not(!) expected.
> I implemented the same methods in a matlab inhouse code & observed perfect 
> convergence rates of 2&2 for v using a direct solver.
>
> The stabilization parameter is in both cases (matlab or deal.ii) chosen as 
> tau = 0.1*h^2/viscosity 
> with 
> double h = std::pow(cell->measure(), 1.0/ static_cast(dim)) 
> or h = cell->diameter() giving similar results.
>
> Could the observed behaviour be caused by applying an iterative solver*? 
> (using ReductionControl and a reduction factor of 1e-13 which is really low)
> * Block-triangular Schur-complement-based approach, similar as in step-55 
> suggested in the further possibilities for extension, basically using 
> S~-(1/viscosity)*Mass_pressure giving 20-40 iterations up to 3mio dofs 
> tested.
>
> I have been checking the code for a week now, and i really doubt, that the 
> integration of just tau*grad(p)*grad(p) is wrong (again, the laplace drops 
> out for a rectangular grid).
> Just to check, i used a manufactured solution from Poisuille flow and 
> added the laplacian from PSPG to the rhs, giving me the exact (linear) 
> solution in the pressure and quadratic convergence in the velocity, thus I 
> assume, that the grad(p)grad(q) term added is correctly implemeted. (exact 
> solution meaning, that i get an L2 error of err<1e-9 for any number of 
> elements used)
>
> Anyone has ever experienced this or has anyone some tipps for further 
> debugging?
> Thanks for taking the time to read the post, any help or comment would be 
> greatly appreciated!
>
> Greetings from Austria,
> Richard
>

-- 
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 

Re: [deal.II] discontinous contour over elements

[Bruno Blais]

Do you use Tecplot to calculate vorticity from the velocity field or do you 
calculate the vorticity from your code, and then visualize it from tecplot?
The way deal.II visualizes gradients (or vorticity in this case) is the 
correct way it should be done, because it is visualized on an "element 
basis", ensuring that the real discontinuities in the fields are seen. Many 
other software will smoothen the discontinuous field by using some sort of 
averaging and then output these values at the nodes. This gives the 
illusion that the fields are continuous, when they are not in reality (the 
averaging introduces the continuity).
Consequently, this really depends on the procedure you are using to 
visualize your other results.


On Saturday, 18 January 2020 19:03:01 UTC-5, David Eaton wrote:
>
> I just use tecplot directly visualize the results. The vorticity contour 
> from my simple code is continuous, and the results from deal.II is 
> discontinuous (without L2 projection). Is it possible that the direct 
> solver in Intel mkl did a similar projection step internally?
>
> --
> *From:* Wolfgang Bangerth >
> *Sent:* Sunday, January 19, 2020, 12:49 AM
> *To:* David Eaton; dea...@googlegroups.com 
> *Subject:* Re: [deal.II] discontinous contour over elements
>
> On 1/18/20 9:25 AM, David Eaton wrote:
> > Thank you for your explanations. Basically I formed a weak form of the 
> PDE for 
> > one element and numerically integrate it at the Gaussian points based on 
> the 
> > interpolation from the local nodes. Subsequently, I assemble the weak 
> forms 
> > from all elements into a global system matrix based on a local-to-global 
> > mapping of the nodes. After applying the boundary conditions, I solve 
> this 
> > linear system using a linear solver in Intel mkl.
>
> Right -- that gives you the solution (u,p) of the problem. But then what 
> do 
> you do to visualize whatever it is that you find is/isn't discontinuous?
>
> 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/95374e76-952b-4b1e-83dc-e784c9c0f472%40googlegroups.com.

Re: [deal.II] discontinous contour over elements

[Bruno Blais]

An easy way to carry the projection that Wolfgang suggested is to use an L2 
projection.
The L2 projection matrix is only a mass matrix and your right hand side is 
constructed by the integral of multiplication of the variable you want to 
project with the test function. Generally, this matrix is very very easy to 
invert. This will yield you the C0 representation of your discontinuous 
field (vorticity) such that the error between your C0 projection and your 
original field is minimized at the position of the gauss points.

This is a procedure we use to set the initial conditions in Lethe when the 
initial condition is a complex function (for instance a Taylor-Green 
vortex).


On Wednesday, 15 January 2020 11:22:52 UTC-5, David Eaton wrote:
>
> I understand the C0 element is piecewise linear across elements. However, 
> I did not experience the same issue in my own C++ code while I use C0 
> element with the Petrov Galerkin stabilization terms. Actually, I am very 
> confused at this point. How could I get rid of it while using C0 element?
>
> Thanks
> D.
>
> --
> *From:* dea...@googlegroups.com   > on behalf of Wolfgang Bangerth  >
> *Sent:* Thursday, January 16, 2020, 12:05 AM
> *To:* dea...@googlegroups.com 
> *Subject:* Re: [deal.II] discontinous contour over elements
>
> On 1/14/20 10:04 PM, David Eaton wrote:
> > 
> > Thank you for your suggestions. I am going to take a look at Lethe and 
> compare 
> > with my implementation. In stabilized formulation, I used quadrilateral 
> > element, instead of P2 P1 Taylor-Hood element. The used element is only 
> C0 
> > element. I also did not expect such a discontinuity between elements. 
> Although 
> > I use P2 P1 Taylor-Hood element without stabilization terms, the 
> discontinuity 
> > is still there. Probably I made mistakes somewhere  in setup. I also 
> suspect 
> > that my solution is not converged. After taking a small relative 
> tolerance 
> > 1e-8, the discontinuity still appears.
>
> David, you did not understand what we were saying: If you use C0 elements 
> (think, piecewise linear) and you take derivatives to compute the 
> vorticity, 
> then you automatically get a discontinuous function. That has nothing to 
> do 
> with stabilization, solver tolerances, etc. It's just a consequence of the 
> fact that C0 elements and their shape functions have kinks and 
> consequently 
> their derivatives are discontinuous.
>
> 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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/f03c64fc-3736-340a-1cdc-f24749fb413f%40colostate.edu
> .
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/bbe79d0e-d423-4df8-8f83-e4e934fedd7c%40googlegroups.com.

Re: [deal.II] discontinous contour over elements

[Bruno Blais]

Dear David,
How are you calculating the vorticity?
As Wolfgang and Praveen have mentioned, if you are using the 
DataPostProcessor, then this will use your shape functions to calculate the 
vorticity. However, your P2-P1 elements are only C0 continuous. 
Consequently, your vorticity can possibly be inherently discontinuous at 
the element edges. However, I am surprised that you obtain such strong 
discontinuity. In our code based on deal.ii (
https://github.com/lethe-cfd/lethe) we implement stabilized formulations 
for the NS equations and the vorticity results for such cases are very 
smooth (even when represented using discontinuous shape functions.

Have you established the convergence of your code using manufactured 
solution? This is the first thing I would suggest. You can look at the 
applications_tests of lethe for some examples of easy manufactured 
solutions for the Incompressible Navier-Stokes equations. There are also 
common books that treat this issue (for instance : 
https://www.amazon.ca/Verification-Validation-Scientific-Computing-Oberkampf/dp/0521113601
)

Please feel free to reach out if you have more questions.
Best
Bruno




On Monday, 13 January 2020 20:29:53 UTC-5, David Eaton wrote:
>
> Thanks Wolfgang and Praveen for providing suggestions. I have tried to 
> debugging this code for a while. I have attached this simple code on this 
> email. I followed the instructions in tutorial closely. Hopefully, anyone 
> could give some suggestions.
>
> Best
> D. 
> --
> *From:* dea...@googlegroups.com   > on behalf of Wolfgang Bangerth  >
> *Sent:* Tuesday, January 14, 2020 6:24 AM
> *To:* dea...@googlegroups.com   >
> *Subject:* Re: [deal.II] discontinous contour over elements 
>  
> On 1/12/20 9:17 PM, David Eaton wrote:
> > My inflow condition is uniform. This formulation and mesh is tested in a 
> > simple C++ code without library. The  large mesh near the inflow does 
> not give 
> > this problem.
> > Yes. I am using C0 element. I did calculation using tecplot. However, 
> the 
> > results from a my C++ code does not give this issue either. Just now, I 
> check 
> > the formulation again. Although I use Q2Q1 Taylor-Hood element without 
> any 
> > stabilization, these issues are still happening.
>
> David -- we don't really know what formulation you are using, how you are 
> implementing it, what you are comparing against, and a number of other 
> factors.
>
> If you have a formulation that computes u,p, and you are plotting the 
> vorticity, you need to expect that the isocontours are discontinuous for 
> the 
> reasons Praveen already stated. If you are getting results that make no 
> sense 
> to you, then the first step would be to ensure that your program is 
> converging 
> as expected. To do this, choose a solution that you know and compute the 
> error 
> norm; then ensure that the program yields error norms that decrease as 
> expected with mesh refinement.
>
> 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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/e1b55e95-a139-d62e-e786-61393cf84dc7%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5843231a-5b39-4f47-8b8e-ca3061959295%40googlegroups.com.

[deal.II] Re: PhD Position (deal.II related) in Computational Fluid Dynamics, Polytechnique Montreal, Montreal, Canada

[Bruno Blais]

Dear Arun,
The project position is filled, but you can always apply spontaneously by 
sending me an email.
Best
Bruno


On Wednesday, 11 December 2019 01:47:23 UTC-5, arun gk wrote:
>
> hi sir.. 
> I would like to apply for this position.
> i would like to know if the vacancy is still there..!?
>
>
> On Saturday, 17 August 2019 00:58:58 UTC+5:30, Bruno Blais wrote:
>>
>> Dear All,
>> We currently have an open position in my research group regarding the 
>> computational modeling of microwave depolymerization reactors using 
>> deal.II. These reactors are part of a key technology that enables the 
>> generation of virgin styrene monomers from contaminated polystyrene using 
>> very novel microwave catalysis technique. This enables the circular economy 
>> of polystyrene plastic waste. The position is fully funded and does not 
>> require the student to carry out TA tasks.
>>
>> A full description is available on the university website: 
>>
>>
>> https://www.polymtl.ca/expertises/sites/expertises.amigow.polymtl.ca/files/2019_polymtl_blais_pyrowave_exp_eng.pdf
>>
>>
>> https://www.polymtl.ca/expertises/en/sustainable-production-chemicals-understanding-microwave-reactor-engineering-hydrodynamics
>>
>> Best regards!
>> Bruno
>>
>

-- 
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/5011e118-1241-415d-8523-d68c29455670%40googlegroups.com.

Re: [deal.II] Getting the normal vector of a cell in where dim

[Bruno Blais]

Wolfgang,
Thanks, this is exactly what I needed at exactly the right place...
I feel slightly dumb right now not having found it :(
Thanks for the quick answer.
Bruno

On Monday, 18 November 2019 12:27:44 UTC-5, Wolfgang Bangerth wrote:
>
>
> Bruno, 
>
> > I am currently working with triangulations that do not have matching 
> > dim-spacedim. For instance a Triangulation<2,3>. For a specific 
> > calculation, I would need the normal vector of the cell of the 
> > triangulation (so the vector perpendicular to the current cell at the 
> > gauss points). Is there any function that would allow me to calculate 
> > this vector? 
> > I have looked at the cellAcessor class and the GeometryInfo struct, but 
> > so far the majority of the content seems geared towards gathering the 
> > normal of the faces of a cell. In this case, I am interested in the 
> > normal of the cell itself, not its faces. 
>
> Yes, you're looking at the wrong place. I bet you need this at the 
> quadrature points, and there the place to look is the FEValues class. It 
> has a member normal() that does what you need when you use it via an 
> FEValues (rather than FEFaceValues) object: 
>
>
> https://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#ac25ec6835799c3b6c7c842f8acb05eb3
>  
>
> 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/d8b25a63-f34d-4172-ba04-37f6649925ae%40googlegroups.com.

[deal.II] Getting the normal vector of a cell in where dim

[Bruno Blais]

Hello,
I am currently working with triangulations that do not have matching 
dim-spacedim. For instance a Triangulation<2,3>. For a specific 
calculation, I would need the normal vector of the cell of the 
triangulation (so the vector perpendicular to the current cell at the gauss 
points). Is there any function that would allow me to calculate this vector?
I have looked at the cellAcessor class and the GeometryInfo struct, but so 
far the majority of the content seems geared towards gathering the normal 
of the faces of a cell. In this case, I am interested in the normal of the 
cell itself, not its faces.

In case this is not readily accessible, would you have a suggestion on how 
to calculate it in an efficient manner? I am thinking that a cross product 
between two vectors of the cell would do the job, but I am afraid that this 
would not work for non Q1 mapping.
Thanks you very much :)!
Bruno

-- 
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/b0ab35c0-8d15-4cf1-8afb-480d3f3999e3%40googlegroups.com.

[deal.II] Re: Convergence rate

[Bruno Blais]

A quick question, since you are working on a sphere, are you specifying a 
mapping of the same order as your phi?

On Wednesday, 9 October 2019 08:57:45 UTC-4, Félix Bunel wrote:
>
> Hello everyone.
>
> I'm having some trouble to understand the convergence rate i'm observing 
> in my code.
>
> Here is what i'm solving :
>
> - I'm in 2D on a round mesh.
> - I'm solving a simple Poisson equation on this mesh for a variable named 
> Phi the solution is known for this and is 1-x^2-y^2
> - With this solution phi I'm then solving a Stokes equation that has 
> special terms that depends on phi.
>
> For the stokes problem i'm using the usual mixed fe element as such :
> stokes_fe(FE_Q<2>(2), 2,
>   FE_Q<2>(1), 1),
> So second order for the speeds and first order for the pressure (just like 
> in the boussinesq problem from the tutorials.
>
> For the poisson/phi problem, i'm using FE_Q<2>  also
> initialized as such :
> phi_fe (2),
> So second order.
>
> For a special case, I have a known solution which is u=v=0 and p of the 
> form 1-x^2-y^2 (just like phi)
> And i have solved this on multiple refinement cycle which gives me 
> different number of dofs and cellsize.
>
> The thing is, when I plot the error as a function of the maximum cell 
> diameter, I get a quadratic convergence rate for the speed, and a linear 
> rate for P and phi.
> What surprises me is that I don't have a quadratic convergence rate for my 
> poisson/\phi problem even though I used 
> phi_fe (2),
>
> My norm is the L2 norm which is simply sqrt(sum(error**2))). I'm computing 
> it with python after exporting the solution as a gpl file.
> Here are the graphs :
>
>
> I have also tried 
> phi_fe (1),
>
> But in this case I don't even have convergence for the speed.
>
>
>
> In my integration part, I have always used 
> QGauss<2>  quadrature_formula(4);
> which is equal to degree+2.
>
> So here are my questions :
>
> 1- Is this convergence rate the one i'm expecting for a poisson equation 
> (phi problem) ?
>
> 2- If no, any idea why I don't have the correct convergence rate ?
>
> 3- In the tutorial step 31, a degree of 2 is used for the temperature 
> (which is very similar to my phi), is there some reason for that choice 
> (just like using degree+1 for the speed and degree for the pressure in a 
> stokes problem...)
>
> Thanks in advance.
>

-- 
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/bf2f032d-3aea-4149-a0b9-c91ecfe176b4%40googlegroups.com.

[deal.II] Re: About Darcy-Brinkman-Forchheimer equation discretization

[Bruno Blais]

Hello, 

It depends of the value of the Reynolds number and the gradients of K with 
respect to x or t, but generally the last two terms do not generally pose 
problems. The first additional term to the right leads to a mass matrix, 
which is well conditioned.
The second term itself is trickier. If you want to use an analytical 
Jacobian formulation and Newton's method, you will need to calculate the 
Frechet derivative of the velocity magnitude. You can also use a Picard 
iteration for this term, which will greatly simplify the expression of the 
Jacobian at the cost more Newton iteration.
Issues generally arise when you have jumps the value of K in space. Then 
generally it is better to use some sort of upwinding to prevent 
oscillations in the velocity field (i.e SUPG).



On Friday, 27 September 2019 22:41:22 UTC-4, FU wrote:
>
> Hi，
> I want to solve the problem about Darcy-Brinkman-Forchheimer equations, 
> but don't know how to discretizate this equation.
>
> [image: Darcy equation.png]
>
> This equation has a similar N-S equation. But the discretization of the 
> last item of the equation and the programming statements are somewhat 
> unclear.
>
>
> Yours，
>
> FU
>

-- 
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/c5bd3188-e6b8-4338-92f3-3a0ec422bb03%40googlegroups.com.

[deal.II] Re: Cmake-GUI error

[Bruno Blais]

Dear Gary,
I use the cmake curse gui a lot (which I guess is what you mean by cmake 
gui).
I have never had any problem installing dealii using it. It should work.


On Wednesday, 18 September 2019 17:26:12 UTC-4, Gary Roach wrote:
>
>
>
> On Tuesday, September 17, 2019 at 3:13:54 PM UTC-7, GaryR wrote:
>>
>> I keep getting the following error when running cmake-gui.
>>
>>
>> CMake Error at /home/gary/Deal.II/share/deal.II/macros/
>> macro_deal_ii_setup_target.cmake:64 (INCLUDE):
>>   INCLUDE could not find load file:
>>
>>
>> /home/gary/Deal.II/lib/cmake/deal.II/deal.IITargets.cmake
>> Call Stack (most recent call first):
>>   /home/gary/Deal.II/share/deal.II/macros/macro_deal_ii_invoke_autopilot.
>> cmake:55 (DEAL_II_SETUP_TARGET)
>>   CMakeLists.txt:39 (DEAL_II_INVOKE_AUTOPILOT)
>>
>>
>> I've tried mucking my way through the CMakeLists.txt. and the 
>> setup_target_.cmake files. I can find the place where the error occurs 
>>  but 
>> am not familiar enough with the coding to make sense of the problem. 
>>
>> I am running Debian Buster with a KDE Desktop.
>> AMD 4 core processor with 16GB of ram.
>> Using Deal.II-9.1.1
>> Using Cmake-GUI-3.13.2
>> Using Eclipse 2019.6 (4.12.0)
>>
>> I have Cmake-GUI set up to produce an Eclipse compatible build.
>>
>> Any help will be sincerely appreciated.
>>
>> GaryR
>>
>
>
>  Sorry about a couple of things. I mistakenly posted the reply with a 
> couple of errors and tried to delete it and re-post. Obviously it didn't 
> work. The CDT3 was a typo. I'm taking your advise and trashing the whole 
> thing and starting over with no GUI.
>
> Thanks for you timely response
>
> GaryR
>
> On Tuesday, September 17, 2019 at 3:13:54 PM UTC-7, GaryR wrote:
>>
>> I keep getting the following error when running cmake-gui.
>>
>>
>> CMake Error at /home/gary/Deal.II/share/deal.II/macros/
>> macro_deal_ii_setup_target.cmake:64 (INCLUDE):
>>   INCLUDE could not find load file:
>>
>>
>> /home/gary/Deal.II/lib/cmake/deal.II/deal.IITargets.cmake
>> Call Stack (most recent call first):
>>   /home/gary/Deal.II/share/deal.II/macros/macro_deal_ii_invoke_autopilot.
>> cmake:55 (DEAL_II_SETUP_TARGET)
>>   CMakeLists.txt:39 (DEAL_II_INVOKE_AUTOPILOT)
>>
>>
>> I've tried mucking my way through the CMakeLists.txt. and the 
>> setup_target_.cmake files. I can find the place where the error occurs 
>>  but 
>> am not familiar enough with the coding to make sense of the problem. 
>>
>> I am running Debian Buster with a KDE Desktop.
>> AMD 4 core processor with 16GB of ram.
>> Using Deal.II-9.1.1
>> Using Cmake-GUI-3.13.2
>> Using Eclipse 2019.6 (4.12.0)
>>
>> I have Cmake-GUI set up to produce an Eclipse compatible build.
>>
>> Any help will be sincerely appreciated.
>>
>> GaryR
>>
>>

-- 
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/cea45c64-7685-4a20-b668-45325e651dff%40googlegroups.com.

Re: [deal.II] Stabilization for Q1Q0 Stokes by Bochev [2006]

[Bruno Blais]

I second Wolfgang comment on the fact that Q1Q1 is not difficult to 
implement. You can also scale it to arbitrary Qn-Qn elements if you are 
interested in higher order.
We have implemented such an approach in our code based on dealii (
https://github.com/lethe-cfd/lethe)
Q1Q1 is already very diffusive (which is why it is so robust I guess), so I 
am not sure that going with Q1-Q0 to save a few pressure degree of freedom 
is actually worth it.
Best!
Bruno

On Tuesday, 10 September 2019 00:03:27 UTC-4, Wolfgang Bangerth wrote:
>
> On 9/9/19 1:57 AM, Richard Schussnig wrote: 
> > 
> > FINALLY, MY QUESTIONS: 
> > 
> > Using the Q1Q1, I would in the end (FSI) need to come up with a space 
> made 
> >  from Q1 elements with a discontinuity at the interface - which shall be 
> > realized using different material_id(). - how may I do that other than 
> > using a FE_DGQ space for the pressure and enforce continuity 'manually' 
> > through a giant ConstraintMatrix? 
>
> That's inefficient, of course :-) I assume that your interface is in the 
> interior of the domain? In that case, take a look at step-46, where 
> solution 
> variables only live on certain cells, and are discontinuous at the 
> interface 
> between the two parts of the domain. 
>
>
> > Using the Q1Q0, the main problem is data transfer and 'node searching' 
> in 
> > the parallel case - example: the stabilization matrix from cell 16 has 
> > pressure dof 45 and shares edges or maybe only a single vertex (!) with 
> > cells with pressure dofs 1 2 3 4 5. The cell matrix for the projection 
> from 
> > Q0(dc) to Q1(c) is an area-weighted sum of the pressures on the cells 
> > touching the vertex of the support of the matching bilinear function, 
> > therefore we get a 6x6 local matrix and entries into all 'touching' 
> cells. 
>
> Yes, you'd have to create a map that for each vertex gives you a list of 
> all 
> adjacent cells. I think I recall that there is a function in GridTools for 
> this, though. 
>
>
> > Since these cells are not only the direct neighbors of the current cell, 
> > things may get complicated quite fast, if we consider the 3d case with 
> > hanging nodes, but on the other hand side, looping in the element loop 
> over 
> > all elements again(!) to check the vertex_index() is extremely slow. 
>
> Yes, you'd reverse this approach by looping over all vertices first, and 
> then 
> in this loop over all adjacent cells. 
>
>
> > Do you know of any better-fitting stabilizations for the Q1Q0 pair? Or 
> do 
> > you think there are better options around? 
>
> Q1-Q1 is a pretty good method, and not very difficult to implement. I'll 
> note 
> that Q1-Q0 *sounds* like a good idea, but has a very low convergence rate 
> and 
> so will not yield particularly good accuracy if that's what you actually 
> care 
> about. Of course, Q2-Q1 is the standard for good reasons. 
>
> 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/169fb0a3-a042-4a01-b25f-2fd46c20bde9%40googlegroups.com.

Re: [deal.II] Difference between BlockDynamicSparsityPattern and TrilinosWrappers::BlockSparsityPattern

[Bruno Blais]

Dear Wolfgang,
Thank you, everything is clear now and i managed to accomplish what I 
wanted. 
Thanks!
Bruno


On Tuesday, 3 September 2019 19:09:00 UTC-4, Wolfgang Bangerth wrote:
>
>
> Bruno, 
>
> > Is there a different between how DynamicSparsityPatterns and 
> > BlockDynamicSparsityPatterns behave? 
>
> The latter is just an array of the former. Under the hood, every block 
> is simply a DynamicSparsityPattern that can be initialized in the same 
> way one always does. 
>
>
> > When you look at step-40, which is the first "MPI" step, the way the 
> > sparsity pattern is made is 
> > (
> https://dealii.org/current/doxygen/deal.II/step_40.html#LaplaceProblemsetup_system)
>  
>
> > | 
> > DynamicSparsityPattern 
> > <
> https://dealii.org/current/doxygen/deal.II/classDynamicSparsityPattern.html>dsp(locally_relevant_dofs);
>  
>
> > DoFTools::make_sparsity_pattern 
> > <
> https://dealii.org/current/doxygen/deal.II/group__constraints.html#gad93530ee35c780e9ef7bc5e4856039df>(dof_handler,dsp,constraints,false);
>  
>
> > SparsityTools::distribute_sparsity_pattern 
> > <
> https://dealii.org/current/doxygen/deal.II/namespaceSparsityTools.html#ae2c7bdbdb62642f60d60087e4cb6195f>(
>  
>
> > dsp, 
> > dof_handler.n_locally_owned_dofs_per_processor(), 
> > mpi_communicator, 
> > locally_relevant_dofs); 
> > | 
> > 
> > My understanding was that you were only making the sparsity pattern for 
> > your own locally relevants dofs and the distribution step would make 
> > sure that everything was coherent. 
>
> Correct. 
>
>
> > However, the sparsity pattern is made 
> > in a completely different way in step-32 
> > (https://dealii.org/current/doxygen/deal.II/step_32.html) and the 
> > TrilinosWrappers:BlockSparsityPattern is used. 
>
> Correct. In step-40, we use PETSc, which has no sparsity pattern data 
> structure of its own, and so we need to initialize the PETSc matrices 
> with our own. In contrast, Trilinos has its own sparsity pattern 
> classes, and so to initialize a Trilinos matrix, we need to build one of 
> their sparsity patterns. (Or block patterns, as it may be -- which again 
> is just an array of patterns.) 
>
>
> > Consequently, I was a bit confused as to the difference between these 
> > approaches. It seems like the second one is necessary for Block 
> matrices? 
> > I will take a deeper look into this, but I must say my understanding of 
> > that point right now is relatively poor :( 
>
> The difference has nothing to do with blocks and everything with whether 
> you base your linear algebra on PETSc or Trilinos. 
>
> Does this help? 
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d7d6b2fe-d2ca-4ea4-9c21-b19953cf2524%40googlegroups.com.

Re: [deal.II] Difference between BlockDynamicSparsityPattern and TrilinosWrappers::BlockSparsityPattern

[Bruno Blais]

Dear Wolfgang,
Thank you very much for your message.
Is there a different between how DynamicSparsityPatterns and 
BlockDynamicSparsityPatterns behave?
When you look at step-40, which is the first "MPI" step, the way the 
sparsity pattern is made is (
https://dealii.org/current/doxygen/deal.II/step_40.html#LaplaceProblemsetup_system
)
DynamicSparsityPattern 
 
dsp(locally_relevant_dofs);
 DoFTools::make_sparsity_pattern 

(dof_handler, dsp, constraints, false);
 SparsityTools::distribute_sparsity_pattern 

(
   dsp,
   dof_handler.n_locally_owned_dofs_per_processor(),
   mpi_communicator,
   locally_relevant_dofs);

My understanding was that you were only making the sparsity pattern for 
your own locally relevants dofs and the distribution step would make sure 
that everything was coherent. However, the sparsity pattern is made in a 
completely different way in step-32 (
https://dealii.org/current/doxygen/deal.II/step_32.html) and the 
TrilinosWrappers:BlockSparsityPattern is used.

Consequently, I was a bit confused as to the difference between these 
approaches. It seems like the second one is necessary for Block matrices?
I will take a deeper look into this, but I must say my understanding of 
that point right now is relatively poor :(
Thanks!
Bruno




On Friday, 30 August 2019 23:41:05 UTC-4, Wolfgang Bangerth wrote:
>
>
> Bruno, 
> I don't quite recall if we ever used the BlockDynamicSparsityPattern in a 
> parallel context. For sure, the way you're initializing it implies that 
> every 
> process allocates the memory for all DoFs as it's not given the 
> information 
> about locally_relevant_dofs. I'd have to look up whether there is a way to 
> do 
> that. On the other hand, I would assume that the Trilinos version has a 
> way to 
> designate what are the locally relevant dofs on each process. 
>
> 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/3d5e5349-ef9a-4a8e-8922-bcd3256caa08%40googlegroups.com.

[deal.II] Difference between BlockDynamicSparsityPattern and TrilinosWrappers::BlockSparsityPattern

[Bruno Blais]

Hello,
I am currently working on a parallel implementation of step-57, thus I am 
learning to live with BlockVectors, BlockMatrices and BlockSparsityPatterns 
in parallel. 
Originally, I thought that I could make my sparsity pattern the following 
way (i.e as in step-57, but distributing it afterwards) :
std::vector block_component(dim+1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise (dof_handler, block_component);
dofs_per_block.resize (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component
);
unsigned int dof_u = dofs_per_block[0];
unsigned int dof_p = dofs_per_block[1];

locally_owned_dofs.resize(2);
locally_owned_dofs[0] = dof_handler.locally_owned_dofs().get_view(0, dof_u);
locally_owned_dofs[1] =dof_handler.locally_owned_dofs().get_view(dof_u, 
dof_u + dof_p);

IndexSet locally_relevant_dofs_acquisition;
DoFTools::extract_locally_relevant_dofs(dof_handler, 
locally_relevant_dofs_acquisition);
locally_relevant_dofs.resize(2);
locally_relevant_dofs[0] = locally_relevant_dofs_acquisition.get_view(0, 
dof_u);
locally_relevant_dofs[1] = locally_relevant_dofs_acquisition.get_view(dof_u, 
dof_u + dof_p);

... Place where I make my constraints ...

BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_sparsity_pattern(dof_handler, dsp, nonzero_constraints);

SparsityTools::distribute_sparsity_pattern(
 dsp,
 dof_handler.locally_owned_dofs_per_processor(),
 mpi_communicator,
 locally_relevant_dofs_acquisition);

system_matrix.reinit(dsp);
pressure_mass_matrix.reinit(dsp.block(1, 1));

When launching in sequential - Debug I have no issue and my solver works 
perfectly. When launching in parallel, within my pre-conditioner, I get an 
error such as :
The violated condition was: 
in.trilinos_partitioner().SameAs(m.DomainMap()) == true
Additional information: 
Column map of matrix does not fit with vector map!

Clearly, I am doing something wrong with my sparsity pattern since it 
appears the block of my vector and the block of my matrix are incompatible 
in size.
I have found that step 32 uses a TrilinosWrappers::BlockSparsityPattern 
instead of a BlockDynamicSparsityPattern. Is this what I should do with my 
case?
I am unsure what is the distinction between what I am doing right now and 
what the TrilinosWrappers::BlockSparsityPattern would do?

If there is a lack of information, I can post a link to the code which is, 
regretfully, quite big.
Best
Bruno

-- 
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/f1b8455d-5a5e-4623-9b70-b4baf87e2435%40googlegroups.com.

[deal.II] Solving a TrilinosWrappers::BlockSparseMatrix using Trilinos iterative solver (step-57)

[Bruno Blais]

Hello,
I am currently trying to work on a "HPC-ready with Trilinos" version of 
step-57 for steady-state solution of the Navier-Stokes equation since I am 
curious about compairing it to our stabilized solver. I think that it  
could give much better results for steady-state.
The approach detailed in step-57  uses Block matrices instead of regular 
matrices since it preconditions the pressure block differently than the 
momentum block.
>From my understanding, this requires eventually the solution of a 
BlockSparseMatrix. However, looking at the documentation for 
TrilinosWrappers::SolverBase, I see that there are no call to solve that 
can take as an argument a BlockSparseMatrix, they all take a regular 
matrix. Consequently, I am not sure how I should proceed. Is this something 
that is lacking in current implementation or there is a work around that I 
am unaware of?

What I am trying to do initially is to port step-57 
(https://www.dealii.org/current/doxygen/deal.II/step_57.html) to use 
Trilinos for the linear algebra. Consequently, what would be changing from 
step-57 is mostly related to changes in the BlockSchurPreconditioner class 
of step-57. Hope that is clear :).

Thanks!
Bruno

-- 
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/b27c883e-4bc2-4ba6-837d-b58f7754f142%40googlegroups.com.

[deal.II] Re: PhD Position (deal.II related) in Computational Fluid Dynamics, Polytechnique Montreal, Montreal, Canada

[Bruno Blais]

As always, it seems I have posted the wrong link.
The link to the PhD related project is :

https://www.polymtl.ca/expertises/sites/expertises.amigow.polymtl.ca/files/2019_polymtl_blais_pyrowave_num_eng_1.pdf

Best!
Bruno

On Friday, 16 August 2019 15:28:58 UTC-4, Bruno Blais wrote:
>
> Dear All,
> We currently have an open position in my research group regarding the 
> computational modeling of microwave depolymerization reactors using 
> deal.II. These reactors are part of a key technology that enables the 
> generation of virgin styrene monomers from contaminated polystyrene using 
> very novel microwave catalysis technique. This enables the circular economy 
> of polystyrene plastic waste. The position is fully funded and does not 
> require the student to carry out TA tasks.
>
> A full description is available on the university website: 
>
>
> https://www.polymtl.ca/expertises/sites/expertises.amigow.polymtl.ca/files/2019_polymtl_blais_pyrowave_exp_eng.pdf
>
>
> https://www.polymtl.ca/expertises/en/sustainable-production-chemicals-understanding-microwave-reactor-engineering-hydrodynamics
>
> Best regards!
> Bruno
>

-- 
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/a5f65860-2274-4ade-b88a-61cb5f1bb92e%40googlegroups.com.

[deal.II] PhD Position (deal.II related) in Computational Fluid Dynamics, Polytechnique Montreal, Montreal, Canada

[Bruno Blais]

Dear All,
We currently have an open position in my research group regarding the 
computational modeling of microwave depolymerization reactors using 
deal.II. These reactors are part of a key technology that enables the 
generation of virgin styrene monomers from contaminated polystyrene using 
very novel microwave catalysis technique. This enables the circular economy 
of polystyrene plastic waste. The position is fully funded and does not 
require the student to carry out TA tasks.

A full description is available on the university website: 

https://www.polymtl.ca/expertises/sites/expertises.amigow.polymtl.ca/files/2019_polymtl_blais_pyrowave_exp_eng.pdf

https://www.polymtl.ca/expertises/en/sustainable-production-chemicals-understanding-microwave-reactor-engineering-hydrodynamics

Best regards!
Bruno

-- 
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/31d54e76-b177-4adc-b8bb-c8a6a2fbe673%40googlegroups.com.

Re: [deal.II] Testing with CTest : Testing the same executable with more than one ".prm" and ".output" combination

[Bruno Blais]

Dear Matthias,
I follow your repository and it works perfectly.
Thanks, I managed to go down from 11 executables to 3... :)!
Best
Bruno


On Wednesday, 14 August 2019 20:04:03 UTC-4, Matthias Maier wrote:
>
>
> On Wed, Aug 14, 2019, at 18:58 CDT, Bruno Blais  > wrote: 
>
> > Dear Matthias, 
> > If I understand correctly, the only constraint of working this way is 
> that 
> > all tests for a single executable must be grouped in a single folder, 
> and 
> > that a single folder can test a single executable? 
> > That's perfect. I will try that tomorrow. 
>
> Exactly. 
>
> > If I manage to make that work, it is exactly what I was looking for. 
>
> Perfect! 
>
> Best, 
> Matthias 
>

-- 
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/ce0230ba-6148-496f-8362-98b6f1a5bb25%40googlegroups.com.

Re: [deal.II] Testing with CTest : Testing the same executable with more than one ".prm" and ".output" combination

[Bruno Blais]

Hahaha, that is amazing!
So much documentation that sometimes documentation becomes its own enemy :P

Thanks for everything! That will definitely be incredibly helpful (not only 
to me, but to others)!

On Wednesday, 14 August 2019 20:00:00 UTC-4, Matthias Maier wrote:
>
>
> On Wed, Aug 14, 2019, at 18:57 CDT, Matthias Maier  > wrote: 
>
> > On Wed, Aug 14, 2019, at 18:55 CDT, Matthias Maier  > wrote: 
> > 
> >> [...] 
> > 
> > This is actually explained here: 
> > 
> >   https://www.dealii.org/current/users/testsuite.html 
>
> AND I apparently created a demo repository showcasing all of that: 
>
> https://github.com/tamiko/dealii-demo 
>
> o_O 
>

-- 
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/0b06f085-6049-4bde-9b56-f0596ff4cf3f%40googlegroups.com.

Re: [deal.II] Testing with CTest : Testing the same executable with more than one ".prm" and ".output" combination

[Bruno Blais]

Dear Matthias,
If I understand correctly, the only constraint of working this way is that 
all tests for a single executable must be grouped in a single folder, and 
that a single folder can test a single executable?
That's perfect. I will try that tomorrow.
If I manage to make that work, it is exactly what I was looking for.

Thanks!
Best
Bruno


On Wednesday, 14 August 2019 19:55:06 UTC-4, Matthias Maier wrote:
>
>
> On Wed, Aug 14, 2019, at 13:38 CDT, Bruno Blais  > wrote: 
>
> > Hello all, 
> > I am re-implementing some tests on our solvers (following 
> > 
> https://www.dealii.org/developer/developers/testsuite.html#layoutaddtests) 
>
> > Right now, it is really easy to test individual applications with a 
> > combination of executable + output + prm file. 
> > We also test with various number of processors. 
> > 
> > However, I was wondering if there was way to test a single executable 
> with 
> > more than one .prm file and corresponding output files 
> > Ex: test.1.prm test.1.output, test.2.prm, test.2.output where test would 
> be 
> > the executable in question. 
> > 
> > Is that something that is feasible in the current framework? 
>
> For example: 
>
> Having in .../tests 
>
>   test_1.prm 
>   test_1.output 
>   test_2.prm 
>   test_2.output 
>
> with the following tests/CMakeLists.txt file: 
>
>   SET(TEST_TARGET my_executable) 
>   DEAL_II_PICKUP_TESTS() 
>
> will define two tests in which the target "my_executable" is run with 
> the corresponding parameter file. For this two work, you will need to 
> have an executable target somewhere else, for example under src: 
>
> src/CMakeLists.txt: 
>
>   ADD_EXECUTABLE(my_executable ) 
>   DEAL_II_SETUP_TARGET(my_executable) 
>
> If you want to run more than one executable, you will need multiple 
> subdirectories (one for each executable). 
>
> Best, 
> Matthias 
>

-- 
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/055bc90a-a3c3-4ff5-af73-1f146e15114e%40googlegroups.com.

[deal.II] Re: Paraview's TemporalInterpolator with adaptive mesh

[Bruno Blais]

Francis (^^)
If you look at the documentation of the temporal interpolator of paraview ( 
https://www.paraview.org/Wiki/Advanced_Animations) you will see that it 
enforces that the support for your data (the mesh) remains constant 
throughout the time series. Consequently, you cannot use it with adaptative 
mesh refinement. 

I think your solution is the best one there is. It might be a good idea, 
for visualization purposes, to project your solution on a refined mesh 
instead of the original coarse one. This way your projection will take into 
account the additional degrees of freedom due to the mesh refinement and 
you will not use too much information in projecting to a static grid. In 
your current case, all of the additional precision due to adaptative mesh 
refinement is lost for the visualization.

Cheers!


On Friday, 12 July 2019 21:51:49 UTC-4, Francis Giraldeau wrote:
>
> I tried SolutionTransfer class, but because the DoFHandler is set in the 
> constructor and assumed to be refined (or coarsen), we cannot use that 
> solution to interpolate to a completely different mesh. Instead, a solution 
> is to use VectorTools::interpolate_to_different_mesh() to a mesh that stays 
> constant:
>
> Vector interpolated_solution;
>
> interpolated_solution.reinit(dof_handler_export.n_dofs());
>
> VectorTools::interpolate_to_different_mesh(dof_handler, solution, 
> dof_handler_export, interpolated_solution);
>
>
> The result can be interpolated using vtkTemporalInterpolator. I uploaded the 
> result at 60fps of step-26 here: https://www.youtube.com/watch?v=1jUUlwOA2fU
>
>
> Cheers,
>
>

-- 
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/2584f6c5-5a8b-4ba7-bdf6-e55769a765ae%40googlegroups.com.

[deal.II] Testing with CTest : Testing the same executable with more than one ".prm" and ".output" combination

[Bruno Blais]

Hello all,
I am re-implementing some tests on our solvers (following 
https://www.dealii.org/developer/developers/testsuite.html#layoutaddtests)
Right now, it is really easy to test individual applications with a 
combination of executable + output + prm file.
We also test with various number of processors.

However, I was wondering if there was way to test a single executable with 
more than one .prm file and corresponding output files
Ex: test.1.prm test.1.output, test.2.prm, test.2.output where test would be 
the executable in question.

Is that something that is feasible in the current framework?
If not, what would be your advice on the way to proceed? I could make 
identical executables that instantiate the same class, but that increases 
linking time for no reasons.

Best
Bruno

-- 
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/59703029-80e0-446f-b47b-b9b04780a7aa%40googlegroups.com.

[deal.II] Error estimation using a function of your solution vector

[Bruno Blais]

Hello everyone, 
I was wondering what was the simplest way to carry out error estimation 
using a function of your solution vector. For example, to carry out mesh 
adaptation using  a Kelly Error estimator but using the vorticity as the 
input variable instead of the velocity vector.
Is there a way to do it in a similar fashion as the way we can post-process 
with DataPostprocessorScalar or DataPostprocessorVector? (which 
is what we use for the vorticity and the q-criterion post-processing)

I have tried looking through the documentation of the Kelly error estimator 
or the steps, but I have not found a similar example.

If such postprocessor cannot be used, what would you suggest to be the best 
way?
I am afraid my knowledge on error estimator is really limited.
Best regards and thank you for your time.
Bruno

-- 
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/495a6a9b-1767-4977-bffc-0c8a7a350870%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Storing additional information at the DOF unrelated to the system you are solving

[Bruno Blais]

Dear Daniel,
I understand, it perfectly makes sense.
I had not thought about using two different FEValues when looping over the 
cells.
Thank you very much!
Best
Bruno


On Monday, July 1, 2019 at 8:52:46 AM UTC-4, Daniel Arndt wrote:
>
> Bruno.
>
> Let's say I solve a first equation for Velocity, and I would like to use 
>> this velocity in another equation for advection-diffusion of say 
>> Temperature.
>> I would first set-up my DOF and FE system for my velocity equation and 
>> solve it and set-up by DOF and FE system for my Temperature.
>> However, how would I be able to set up an additional space to store my 
>> velocity at the nodes that are related to my Temperature FESystem without 
>> introducing additional unknown in my  temperature equations (Since velocity 
>> is known after all) and without assuming that both my temperature and 
>> velocity equation are numbered in the same exact fashion? 
>>
>
> For assemling a bilinear form you typically don't nees the support points 
> of the ansatz functions but their evaluations in the quadrature points.
> In your case, I would simply use two DoFHandler objects (one for each 
> system of equations) and loop over both of them for the temperature 
> equations.
>
> https://www.dealii.org/current/doxygen/deal.II/step_31.html#BoussinesqFlowProblemassemble_temperature_system
>  
> might be helpful.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f14e1cd9-ca06-421f-9b21-b89e2c1c04ea%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Storing additional information at the DOF unrelated to the system you are solving

[Bruno Blais]

Dear Daniel,
I am sorry, my question truly was unclear.
Let's say I solve a first equation for Velocity, and I would like to use 
this velocity in another equation for advection-diffusion of say 
Temperature.
I would first set-up my DOF and FE system for my velocity equation and 
solve it and set-up by DOF and FE system for my Temperature.
However, how would I be able to set up an additional space to store my 
velocity at the nodes that are related to my Temperature FESystem without 
introducing additional unknown in my  temperature equations (Since velocity 
is known after all) and without assuming that both my temperature and 
velocity equation are numbered in the same exact fashion? 

Best
Bruno




On Saturday, June 29, 2019 at 1:54:03 PM UTC-4, Daniel Arndt wrote:
>
> Bruno,
>
> Hello everyone, I hope you are well.
>> I have a quick question which I cannot seem to wrap my hear around. I 
>> think it will sound confused, but anyway.
>> Let's say I am solving for a scalar equation. I would have a 
>> triangulation, then a DoFHandler and a FESystem with a multiplicity of one.
>> I can loop over the cells and assemble my equations.
>>
> If you only have a scalar equation, there is no need to use FESystem but 
> you can do that of course.
>  
>
>> Now, let's say I wanted to have an additional property that would be  a 
>> vector defined at the nodes of my triangulation using the positions of the 
>> nodes.
>> However, this vector would not be used to be interpolated at the gauss 
>> points, but would be use to cut the cells and do other geometric 
>> manipulations.
>> I could define another DoFHandler , and FESystem, but then when I loop 
>> over the cells, there would be nothing to make sure that my DOF on both 
>> equations are numbered in the same way.
>> However, I cannot defined another vector using the FESystem of my scalar 
>> equation, because the multiplicity is not the same.
>> What would be the best way to proceed?
>>
> I am not quite sure to understand what you are asking for. The support 
> points only coincide for FE_Q(1) elements with the nodes of the 
> triangulation.
> Normally, you would associate the entries of a Vector with the finite 
> element vector defined by the degrees of freedom (and the ansatz space), 
> but you can also store different information.
> You are talking about two equations. What is the second equation and how 
> is it related to the first one? Can you give an example how you want to use 
> the Vectors?
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/bf313a02-e47a-432b-8d7b-f3d4be42f8fc%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Storing additional information at the DOF unrelated to the system you are solving

[Bruno Blais]

Hello everyone, I hope you are well.
I have a quick question which I cannot seem to wrap my hear around. I think 
it will sound confused, but anyway.
Let's say I am solving for a scalar equation. I would have a triangulation, 
then a DoFHandler and a FESystem with a multiplicity of one.
I can loop over the cells and assemble my equations.

Now, let's say I wanted to have an additional property that would be  a 
vector defined at the nodes of my triangulation using the positions of the 
nodes.
However, this vector would not be used to be interpolated at the gauss 
points, but would be use to cut the cells and do other geometric 
manipulations.
I could define another DoFHandler , and FESystem, but then when I loop over 
the cells, there would be nothing to make sure that my DOF on both 
equations are numbered in the same way.
However, I cannot defined another vector using the FESystem of my scalar 
equation, because the multiplicity is not the same.
What would be the best way to proceed?

TLDR: I would like to store additional information at the position of the 
DOFs, without including them in my vector of unknown, but I would still 
want to be able to know their position.
Is this doable?
Thanks!
Bruno


-- 
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/8b4d0f2f-8889-4d97-bc3d-9d04c93c0167%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Sliding mesh-like approach using dealii

[Bruno Blais]

Dear Earl and Wolfgang,
Thanks for the very thorough answer.
I will definitely look into it this year and try to come up with a 
solution. I think there are many ways to tackle this type of approach using 
the existing tools :)
I'll keep you informed should I find a way to make a decent implementation.
Best!
Bruno


On Sunday, 16 June 2019 21:13:32 UTC-4, Wolfgang Bangerth wrote:
>
>
> Bruno, 
>
> > As some of you might know, sliding mesh approaches are generally used in 
> CFD 
> > simulation of rotating geometries without axial symmetry (for instance, 
> an 
> > impeller with baffles). 
> > This is generally achieved by having two triangulation, one that is 
> rotating 
> > and one that is static. At the interface between the meshes, constraints 
> are 
> > used to "bridge" the two meshes together. Although this generally 
> induces 
> > additional interpolation error, this is generally one of the best way to 
> deal 
> > with turbomachinery. 
> > 
> > Since dealii is already equipped to deal with contact problem, has 
> anybody 
> > ever investigated if a sliding-mesh type of simulation could be carried 
> out 
> > using dealii? 
>
> I don't know, but then there are 1,200+ publications at 
>https://dealii.org/publications.html 
> that might contain something you're looking for. There's now even a search 
> function :-) 
>
> There are going to be two challenges: 
> * Mesh generation. Earl Fairall already commented on that. 
> * Generation of the constraints. That's going to be difficult because 
>you'll have to find quadrature points on the faces of one mesh on the 
>faces of the other mesh, and these sorts of point search algorithms 
>are always expensive (though easy to parallelize). You need this kind 
>of mapping if you want to project between the two mesh surfaces; for 
>interpolation, you'll need to find the location of the support points 
>of the faces of one mesh on the faces of the other mesh, which is the 
>same kind of operation. 
>
> You could consider a mortar approach in which you would have a third mesh 
> at 
> the interface. In that case, you could choose at least one of the meshes 
> involved to be somewhat structured, but in the end, it's probably going to 
> be 
> about as expensive as otherwise. 
>
> I don't know anyone who has implemented the exact kind of application you 
> have, but you might want to look up some of the work done on 
> fluid-structure 
> interaction in deal.II (e.g., by Thomas Wick). I *believe* that they too 
> have 
> to interpolate between different kinds of meshes. There is also the 
> 'nonmatching' namespace that was added by Luca Heltai and coworkers in the 
> last release that helps you deal with overlay meshes -- which would also 
> be a 
> way to do what you're looking at -- the rotating geometry would simply be 
> an 
> overlay to a background mesh. I believe there is also a tutorial program 
> for that. 
>
> Other than that, I have no real pointers. But it's an interesting topic, 
> and 
> if you find ways to implement what you are looking for, please feel free 
> to 
> post solutions here (and/or make small test programs available as code 
> gallery 
> or tutorial program! 
>
> 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/3fdbf6a1-cd04-4e0d-82bd-859a032091e3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Sliding mesh-like approach using dealii

[Bruno Blais]

Hello everyone,
As some of you might know, sliding mesh approaches are generally used in 
CFD simulation of rotating geometries without axial symmetry (for instance, 
an impeller with baffles).
This is generally achieved by having two triangulation, one that is 
rotating and one that is static. At the interface between the meshes, 
constraints are used to "bridge" the two meshes together. Although this 
generally induces additional interpolation error, this is generally one of 
the best way to deal with turbomachinery.

Since dealii is already equipped to deal with contact problem, has anybody 
ever investigated if a sliding-mesh type of simulation could be carried out 
using dealii?

Thanks
Bruno

-- 
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/2c2492a4-cbbb-4c83-b550-f398a0f08358%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Periodic boundary conditions and the AffineConstraints class

[Bruno Blais]

Indeed it was my mistake, sorry :)


On Monday, 6 May 2019 01:48:23 UTC-4, Matthias Maier wrote:
>
>
> On Sun, May  5, 2019, at 22:47 CDT, Bruno Blais  > wrote: 
>
> > I am trying to implement periodic boundary conditions in my CFD solver. 
> > When looking at the documentation for the make_periodicity_constraints I 
> > see that the contraints have to be of type ConstraintsMatrix instead of 
> > type AffineConstraints (or another type) 
>
> > Is this a current issue or it is better to stick with ConstraintsMatrix 
> for 
> > the time being? 
>
> You have probably looked at different release versions for both 
> classes. We have renamed ConstraintMatrix to AffineConstraints (and made 
> ConstraintMatrix a type alias for AffineConstraints). 
>
> The function make_periodicity_constraints will thus of course work with 
> AffineConstraints. 
>
> Best, 
> Matthias 
>

-- 
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/f45306aa-8395-4f5b-aabc-a35e9d22e8f9%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Periodic boundary conditions and the AffineConstraints class

[Bruno Blais]

Hello everyone,
I am trying to implement periodic boundary conditions in my CFD solver.
When looking at the documentation for the make_periodicity_constraints I 
see that the contraints have to be of type ConstraintsMatrix instead of 
type AffineConstraints (or another type)
Is this a current issue or it is better to stick with ConstraintsMatrix for 
the time being?
Thank you for everything,
Bruno


-- 
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] Setting an initial condition by fixing the DOF values themselves

[Bruno Blais]

Dear Wolfgang, this was easier than I thought.
I followed your suggestion and used VectorTools::interpolate.
Everything works as expected.
Thank you very much

Final solution:

 // initialConditionParameters.uvw is the parsed function
 // newton_update is a local vector

  const FEValuesExtractors::Vector velocities (0);

  const FEValuesExtractors::Scalar pressure (dim);

  const MappingQ  mapping (degreeVelocity_,femParameters.qmapping_all);


  VectorTools::interpolate(mapping,

   dof_handler,

   initialConditionParameters.uvw,

   newton_update,

   fe.component_mask(velocities));


Thanks again!


On Thursday, 25 April 2019 09:32:09 UTC-4, Wolfgang Bangerth wrote:
>
> On 4/25/19 6:58 AM, Bruno Blais wrote: 
> > I would like to be able to set-up an initial condition by fixing the 
> value of 
> > the DOF themselves by using their x,y(,z) position with a ParsedFuction. 
> > 
> > Generally, I set the initial condition by using an L2 projection of the 
> > function (for complex function and higher order element I know this is 
> more 
> > appropriate). 
> > However, for some stuff I would like to be able to fix the value of the 
> DOF 
> > directly. 
> > My issue is, through the doxygen, I have not found a way to loop through 
> the 
> > DOF and to extract their position, even more when I have high order 
> element or 
> > when I have a multiplicity (say velocity-vector + pressure). 
> > 
> > What I would be looking into doing would be something similar to : 
> > 
> > loop over local DOF; 
> >get DOF position 
> >get DOF component (u[0], u[1], u[2], p) 
> >calculate parsed function using the point 
> >set the value of the DOF 
> > end 
> > 
> > Is it something that is possible? This is a weird way of iterating since 
> I 
> > know we always loop over the cells. Maybe there is a way to achieve this 
> by 
> > looping over the cells and then looping over the DOF of the cells? This 
> would 
> > be slightly redundant, but I do not care since I am doing this only 
> once. 
>
> You can get the locations of DoFs using 
> DoFTools::map_dofs_to_support_points(). 
>
> But maybe the easier way is to use 
> VectorTools::interpolate_boundary_values() 
> or project() and just pass it a function object that describes the 
> function 
> you want to set in terms of (x,y,z). This could be your ParsedFunction. 
> interpolate_b_v() will then call your parsed function with a bunch of 
> points 
> that happen to correspond to the DoF locations, but it makes the 
> connection 
> internally and you no longer need to worry about which DoF sits where. 
>
> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Setting an initial condition by fixing the DOF values themselves

[Bruno Blais]

Dear Wolfgang,
This was easier than I thought. I followed your suggestion and used 
VectorTools::interpolate
Final solution to set the values of the DOF using a parsed function:

// initialConditionParameters.uvw is my parsed function
// newtonUpdate is my local vector that I use to alter the DOF values


  const FEValuesExtractors::Vector velocities (0);

  const FEValuesExtractors::Scalar pressure (dim);

  const MappingQ  mapping (degreeVelocity_,femParameters.qmapping_all);

  VectorTools::interpolate(mapping, 

   dof_handler,

   initialConditionParameters.uvw,

   newton_update,

   fe.component_mask(velocities));


Thank you very much for the help!

On Thursday, 25 April 2019 09:32:09 UTC-4, Wolfgang Bangerth wrote:
>
> On 4/25/19 6:58 AM, Bruno Blais wrote: 
> > I would like to be able to set-up an initial condition by fixing the 
> value of 
> > the DOF themselves by using their x,y(,z) position with a ParsedFuction. 
> > 
> > Generally, I set the initial condition by using an L2 projection of the 
> > function (for complex function and higher order element I know this is 
> more 
> > appropriate). 
> > However, for some stuff I would like to be able to fix the value of the 
> DOF 
> > directly. 
> > My issue is, through the doxygen, I have not found a way to loop through 
> the 
> > DOF and to extract their position, even more when I have high order 
> element or 
> > when I have a multiplicity (say velocity-vector + pressure). 
> > 
> > What I would be looking into doing would be something similar to : 
> > 
> > loop over local DOF; 
> >get DOF position 
> >get DOF component (u[0], u[1], u[2], p) 
> >calculate parsed function using the point 
> >set the value of the DOF 
> > end 
> > 
> > Is it something that is possible? This is a weird way of iterating since 
> I 
> > know we always loop over the cells. Maybe there is a way to achieve this 
> by 
> > looping over the cells and then looping over the DOF of the cells? This 
> would 
> > be slightly redundant, but I do not care since I am doing this only 
> once. 
>
> You can get the locations of DoFs using 
> DoFTools::map_dofs_to_support_points(). 
>
> But maybe the easier way is to use 
> VectorTools::interpolate_boundary_values() 
> or project() and just pass it a function object that describes the 
> function 
> you want to set in terms of (x,y,z). This could be your ParsedFunction. 
> interpolate_b_v() will then call your parsed function with a bunch of 
> points 
> that happen to correspond to the DoF locations, but it makes the 
> connection 
> internally and you no longer need to worry about which DoF sits where. 
>
> 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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Setting an initial condition by fixing the DOF values themselves

[Bruno Blais]

Hello everyone,
I believe this is a relatively dumb question, but I seem to struggle with 
this basic concept.
I would like to be able to set-up an initial condition by fixing the value 
of the DOF themselves by using their x,y(,z) position with a ParsedFuction.

Generally, I set the initial condition by using an L2 projection of the 
function (for complex function and higher order element I know this is more 
appropriate).
However, for some stuff I would like to be able to fix the value of the DOF 
directly.
My issue is, through the doxygen, I have not found a way to loop through 
the DOF and to extract their position, even more when I have high order 
element or when I have a multiplicity (say velocity-vector + pressure).

What I would be looking into doing would be something similar to :

loop over local DOF;
  get DOF position
  get DOF component (u[0], u[1], u[2], p)
  calculate parsed function using the point
  set the value of the DOF
end

Is it something that is possible? This is a weird way of iterating since I 
know we always loop over the cells. Maybe there is a way to achieve this by 
looping over the cells and then looping over the DOF of the cells? This 
would be slightly redundant, but I do not care since I am doing this only 
once.

Thank you for everything!
Bruno


-- 
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] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

Dear Jean-Paul,
Thank you, everything works. Right now for simple test cases I have managed 
to get interesting results:

   - Going from ILU to AMG leads to speed-up of 4X on the same number of 
   processor.
   -  More importantly, my number of GMRES iteration remains extremely low 
   (below 15) and is practically independent of the number of cells (10 
   iterations with 128x128, 12 with 256x256, etc.). If my multigrid solver 
   classes are not too far behind me, this should be the expected behavior

I need to test it more on more complex test cases, but once this is done, I 
will open a pull request and try to make an interesting tutorial out of it 
this summer.
I just wanted to say thanks for your amazing contribution, these small 
modifications really open up a lot of possibilities in terms of exploring 
the AMG.

With this, I have weeks of fun ahead :)! 

Here is a small final example for the setting of the ILU smoother and 
coarse parameters.

TrilinosWrappers::PreconditionAMG preconditioner;

  std::vector > constant_modes;
  std::vector  velocity_components (dim+1,true);
  DoFTools::extract_constant_modes (dof_handler, velocity_components,
constant_modes);
  TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
  amg_data.constant_modes = constant_modes;

  const boolelliptic=false;
  const boolhigher_order_elements = false;
  const unsigned int  n_cycles = 2;
  const boolw_cycle = true;
  const double  aggregation_threshold = 1e-10;
  const unsigned int  smoother_sweeps = 2;
  const unsigned int  smoother_overlap = 1;
  const booloutput_details = false;
  const char *  smoother_type = "ILU";
  const char *  coarse_type = "ILU";
  TrilinosWrappers::PreconditionAMG::AdditionalData preconditionerOptions(
elliptic,
higher_order_elements,
n_cycles,
w_cycle,
aggregation_threshold,
constant_modes,
smoother_sweeps,
smoother_overlap,
output_details,
smoother_type,
coarse_type
);

  Teuchos::ParameterListparameter_ml;
  std::unique_ptr< Epetra_MultiVector > distributed_constant_modes;
  preconditionerOptions.set_parameters(parameter_ml, 
distributed_constant_modes, system_matrix);
  const double ilu_fill=linearSolverParameters.ilut_fill;
  const double ilu_atol=linearSolverParameters.ilut_atol ;
  const double ilu_rtol=linearSolverParameters.ilut_rtol;
  parameter_ml.set("smoother: ifpack level-of-fill",ilu_fill);
  parameter_ml.set("smoother: ifpack absolute threshold",ilu_atol);
  parameter_ml.set("coarse: ifpack level-of-fill",ilu_fill);
  parameter_ml.set("coarse: ifpack absolute threshold",ilu_atol);
  preconditioner.initialize(system_matrix,parameter_ml);



On Tuesday, 2 April 2019 05:55:26 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> These are the basic steps that you need to take:
>
> 1. Create and initialise a PreconditionAMG::AdditionalData object 
> (“additional_data”) as per usual. Here you set up the constant modes as is 
> needed using DoFTools::extract_constant_modes().
> 2. Create a Teuchos::ParameterList (“parameter_ml”) and std::unique_ptr< 
> Epetra_MultiVector > (“distributed_constant_modes”). Neither of these need 
> be initialised with any data.
> 3. Now call additional_data.set_parameters(parameter_ml, 
> distributed_constant_modes, matrix); .  This initialises “ parameter_ml” 
> and “”distributed_constant_modes” for you, and transcribes the settings 
> stored in “additional_data” into the equivalent “parameters_ml”. The 
>  “distributed_constant_modes” is internally referenced by “ parameters_ml”, 
> so thats why you need that vector lying around.
> 4. Lastly, you initialise the preconditioner using ” parameter_ml” for the 
> “matrix” that you sent into “ additional_data.set_parameters()”, i.e. 
> something like this: const PreconditionAMG prec (matrix, parameter_ml); .
>
> As for the second point, you should be able to find a link to the ML 
> documentation in the description to this function:
>
> https://dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1PreconditionAMG.html#a9703cc100be147a47358cf99fc53beb0
>
> You can also search for this technical report which should outline most if 
> not all of the possible settings:
>
> *TechReport (Gee2006a)*
> Gee, M. W.; Siefert, C. M.; Hu, J. J.; Tuminaro, R. S. & Sala, M. G.
> ML 5.0 Smoothed Aggregation User's Guide 
> *Sandia National Laboratories, * *Sandia National Laboratories, * *2006* 
>
> Best,
> Jean-Paul
>
> On 01 Apr 2019, at 17:53, Bruno Blais > 
> wrote:
>
> Dear Jean-Paul,
> I've compiled the new DEALII version with your changes, but I am unsure of 
> how to use the cha

Re: [deal.II] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

Dear Jean-Paul,
Thank you very much. I will look into this and the documentation.
I will post you the final results.
I greatly appreciate your support and these very nice changes to the 
wrapper :)
Best
Bruno


On Tuesday, 2 April 2019 05:55:26 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> These are the basic steps that you need to take:
>
> 1. Create and initialise a PreconditionAMG::AdditionalData object 
> (“additional_data”) as per usual. Here you set up the constant modes as is 
> needed using DoFTools::extract_constant_modes().
> 2. Create a Teuchos::ParameterList (“parameter_ml”) and std::unique_ptr< 
> Epetra_MultiVector > (“distributed_constant_modes”). Neither of these need 
> be initialised with any data.
> 3. Now call additional_data.set_parameters(parameter_ml, 
> distributed_constant_modes, matrix); .  This initialises “ parameter_ml” 
> and “”distributed_constant_modes” for you, and transcribes the settings 
> stored in “additional_data” into the equivalent “parameters_ml”. The 
>  “distributed_constant_modes” is internally referenced by “ parameters_ml”, 
> so thats why you need that vector lying around.
> 4. Lastly, you initialise the preconditioner using ” parameter_ml” for the 
> “matrix” that you sent into “ additional_data.set_parameters()”, i.e. 
> something like this: const PreconditionAMG prec (matrix, parameter_ml); .
>
> As for the second point, you should be able to find a link to the ML 
> documentation in the description to this function:
>
> https://dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1PreconditionAMG.html#a9703cc100be147a47358cf99fc53beb0
>
> You can also search for this technical report which should outline most if 
> not all of the possible settings:
>
> *TechReport (Gee2006a)*
> Gee, M. W.; Siefert, C. M.; Hu, J. J.; Tuminaro, R. S. & Sala, M. G.
> ML 5.0 Smoothed Aggregation User's Guide 
> *Sandia National Laboratories, * *Sandia National Laboratories, * *2006* 
>
> Best,
> Jean-Paul
>
> On 01 Apr 2019, at 17:53, Bruno Blais > 
> wrote:
>
> Dear Jean-Paul,
> I've compiled the new DEALII version with your changes, but I am unsure of 
> how to use the changes correctly.
> If I understand correctly, you can now set Teuchos parameters after you 
> have initialized the TrilinosWrappers::PreconditionAMG.
> This is achieved by creating a Teuchos::ParameterList object and using it 
> in the set_parameters public function of the 
> Preconditioner::AdditionalParameters.
> So for instance in my case it would be something like:
>
>
> // Still need to figure out the syntax for the additional_parameters
>  Teuchos::ParameterListadditional_parameter("smoother: params"
> );
>  
>
> And then I would set the parameters in the preconditioner additional 
> parameters
>
> I have two things i am unsure:
> 1) For the set the additional parameters I need a 
> std::unique_ptr _constant_modes . I presume 
> this vector is created in the initialization from the std::vector of 
> constant boolean modes? Is there a way to grab it again and reuse it? I do 
> not understand why I need to specify the constant modes for every change of 
> parameters.
> 2) Is there any documentation on the various parameter list. What I would 
> like to change is the atol of the ILU smoother and ILU preconditions. I 
> think this will be in the smoother: params, but I am unsure how to proceed 
> from there. I have found the MULUE documentation, but it seems to be more 
> tailored around their XML interface?
>
>
> Thank you so much,
> Bruno
>
>
>
> On Tuesday, 12 March 2019 05:06:03 UTC-4, Jean-Paul Pelteret wrote:
>>
>> Dear Bruno,
>>
>> I believe that my proposed changes are working now, so you could always 
>> pull that branch from my repository and use it immediately. That might be 
>> the most convenient way for you to move forward.
>>
>> Thanks for being willing to share your findings wrt. AMG parameters. I 
>> think that would be really useful. I also think that encapsulating a 
>> parameter study for NS + AMG/GMG within a tutorial could be interesting and 
>> valuable. Maybe you could open up an issue on the GitHub repository and we 
>> could all discuss the specifics there.
>>
>> Best,
>> Jean-Paul
>>
>> On 11 Mar 2019, at 14:25, Bruno Blais  wrote:
>>
>>  Dear Jean-Paul and Wolfgang,
>> Is it better if I try to go with the Teuchos way or should I wait until 
>> Jean-Paul finishes the merge on his branch ? I think the latter might be a 
>> good option.
>> In all cases, I would be more than glad to share my experience in 
>> parameters tuning for the AMG solver with ILU t

Re: [deal.II] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

Dear Jean-Paul,
I've compiled the new DEALII version with your changes, but I am unsure of 
how to use the changes correctly.
If I understand correctly, you can now set Teuchos parameters after you 
have initialized the TrilinosWrappers::PreconditionAMG.
This is achieved by creating a Teuchos::ParameterList object and using it 
in the set_parameters public function of the 
Preconditioner::AdditionalParameters.
So for instance in my case it would be something like:


// Still need to figure out the syntax for the additional_parameters
 Teuchos::ParameterListadditional_parameter("smoother: params");
 

And then I would set the parameters in the preconditioner additional 
parameters

I have two things i am unsure:
1) For the set the additional parameters I need a 
std::unique_ptr _constant_modes . I presume 
this vector is created in the initialization from the std::vector of 
constant boolean modes? Is there a way to grab it again and reuse it? I do 
not understand why I need to specify the constant modes for every change of 
parameters.
2) Is there any documentation on the various parameter list. What I would 
like to change is the atol of the ILU smoother and ILU preconditions. I 
think this will be in the smoother: params, but I am unsure how to proceed 
from there. I have found the MULUE documentation, but it seems to be more 
tailored around their XML interface?


Thank you so much,
Bruno



On Tuesday, 12 March 2019 05:06:03 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> I believe that my proposed changes are working now, so you could always 
> pull that branch from my repository and use it immediately. That might be 
> the most convenient way for you to move forward.
>
> Thanks for being willing to share your findings wrt. AMG parameters. I 
> think that would be really useful. I also think that encapsulating a 
> parameter study for NS + AMG/GMG within a tutorial could be interesting and 
> valuable. Maybe you could open up an issue on the GitHub repository and we 
> could all discuss the specifics there.
>
> Best,
> Jean-Paul
>
> On 11 Mar 2019, at 14:25, Bruno Blais > 
> wrote:
>
>  Dear Jean-Paul and Wolfgang,
> Is it better if I try to go with the Teuchos way or should I wait until 
> Jean-Paul finishes the merge on his branch ? I think the latter might be a 
> good option.
> In all cases, I would be more than glad to share my experience in 
> parameters tuning for the AMG solver with ILU type of smoother and coarse 
> solve.
>
> On a slightly related note. Hopefully I will have finished the basis of my 
> GLS stabilized navier-stokes solver by June or so. It already works pretty 
> well in parallel with GMRES + ILU, but I want to test more what can be done 
> in terms of multigrid and AMG.
> In a way, it is very reminescent of a similar work done in MOOSE, but 
> using Trilinos, P4est and deal.II. The MOOSE paper is here:
> https://www.sciencedirect.com/science/article/pii/S0965997817310591
>
> I am more than willing to share everything about the solver. Do you 
> believe this could be something that could make an interesting Tutorial? I 
> could wrap-it around the study of the flow around a cylinder or a backward 
> facing step to have some validation. Verification would be done using MMS.
> Best
> Bruno
>
> -- 
> 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 .
> For more options, visit https://groups.google.com/d/optout.
>
>
>

-- 
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] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

Dear Jean-Paul.
That's good, I will test the AMG with the stabilized approach and see what 
I can get.
If the results are interesting (well they already are right now for ILU + 
GMRES) i'll open up an issue on the github for the dicussion on that more 
specialized topic.
Thank you for everything
Bruno

On Tuesday, 12 March 2019 05:06:03 UTC-4, Jean-Paul Pelteret wrote:
>
> Dear Bruno,
>
> I believe that my proposed changes are working now, so you could always 
> pull that branch from my repository and use it immediately. That might be 
> the most convenient way for you to move forward.
>
> Thanks for being willing to share your findings wrt. AMG parameters. I 
> think that would be really useful. I also think that encapsulating a 
> parameter study for NS + AMG/GMG within a tutorial could be interesting and 
> valuable. Maybe you could open up an issue on the GitHub repository and we 
> could all discuss the specifics there.
>
> Best,
> Jean-Paul
>
> On 11 Mar 2019, at 14:25, Bruno Blais > 
> wrote:
>
>  Dear Jean-Paul and Wolfgang,
> Is it better if I try to go with the Teuchos way or should I wait until 
> Jean-Paul finishes the merge on his branch ? I think the latter might be a 
> good option.
> In all cases, I would be more than glad to share my experience in 
> parameters tuning for the AMG solver with ILU type of smoother and coarse 
> solve.
>
> On a slightly related note. Hopefully I will have finished the basis of my 
> GLS stabilized navier-stokes solver by June or so. It already works pretty 
> well in parallel with GMRES + ILU, but I want to test more what can be done 
> in terms of multigrid and AMG.
> In a way, it is very reminescent of a similar work done in MOOSE, but 
> using Trilinos, P4est and deal.II. The MOOSE paper is here:
> https://www.sciencedirect.com/science/article/pii/S0965997817310591
>
> I am more than willing to share everything about the solver. Do you 
> believe this could be something that could make an interesting Tutorial? I 
> could wrap-it around the study of the flow around a cylinder or a backward 
> facing step to have some validation. Verification would be done using MMS.
> Best
> Bruno
>
> -- 
> 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 .
> For more options, visit https://groups.google.com/d/optout.
>
>
>

-- 
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] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

 Dear Jean-Paul and Wolfgang,
Is it better if I try to go with the Teuchos way or should I wait until 
Jean-Paul finishes the merge on his branch ? I think the latter might be a 
good option.
In all cases, I would be more than glad to share my experience in 
parameters tuning for the AMG solver with ILU type of smoother and coarse 
solve.

On a slightly related note. Hopefully I will have finished the basis of my 
GLS stabilized navier-stokes solver by June or so. It already works pretty 
well in parallel with GMRES + ILU, but I want to test more what can be done 
in terms of multigrid and AMG.
In a way, it is very reminescent of a similar work done in MOOSE, but using 
Trilinos, P4est and deal.II. The MOOSE paper is here:
https://www.sciencedirect.com/science/article/pii/S0965997817310591

I am more than willing to share everything about the solver. Do you believe 
this could be something that could make an interesting Tutorial? I could 
wrap-it around the study of the flow around a cylinder or a backward facing 
step to have some validation. Verification would be done using MMS.
Best
Bruno

-- 
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] Controllins fill-in, abs. tol and rel. tol or smoother and coarse solver of TrilinosWrappers::PreconditionAMG

[Bruno Blais]

Hello everyone,
Just a small question.
Looking at the documentation and at the .h of the  
TrilinosWrappers::PreconditionAMG  I have not found a way to access the 
smoother and coarse solver to play with their parameters. Say if you use an 
ILU or ILUT smoother.
When you use the verbose output of the solver you get status report such as:

Smoother (level 2) : IFPACK, type=`ILU',
Smoother (level 2) : post,overlap=1
Smoother (level 2) : level-of-fill=0,rel. threshold=1,abs. threshold=0
Smoother (level 2) : Setup time : 0.00463979 (s)

Coarse solve (level 3) : IFPACK, type=`ILU',
Coarse solve (level 3) : post,overlap=1
Coarse solve (level 3) : level-of-fill=0,rel. threshold=1,abs. threshold=0


Clearly, parameters like fill-level, rel. treshold and abs. threshold could 
be modified. Is there any way to access them and modify them? I have not 
been able to as of now.

I would prefer to keep using the Wrappers as I find their usage to be 
relatively easy and straightforward.
Thank you for everything,
Best
Bruno

-- 
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] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

On Thursday, 7 March 2019 23:35:15 UTC-5, Wolfgang Bangerth wrote:
>
> On 3/7/19 7:45 AM, Bruno Blais wrote: 
> > 
> > However, my GMRES stops very quickly after a certain number of newton 
> > iteration with the following : AztecOO::Iterate error code -4: GMRES 
> > Hessenberg ill-conditioned 
>
> This is surprising and suggests that something is wrong with the matrix, 
> not 
> the right hand side. At least that's what I think this probably means. Did 
> you 
> have the same error when you did not modify the right hand side? What 
> happens 
> if you just artificially create a right hand side that you know is in the 
> range, for example by picking a random vector and multiplying it by the 
> matrix? 
>

I did not have the same error when I did not modify the RHS. Which I also 
find very surprising.
I also found that I did not get the same error if I modify the RHS, but 
that I also modify the Newton correction (if I project my newton correction 
such that corr = corr - P*corr)
I will look at what you suggested, that's an extremely good idea, and get 
back at you.
Thanks!
 

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

Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

Dear Wolfgang,
Sorry for an additional question.
I have tried to implement by myself something similar to what is done in 
ASPECT (I guess?, it is a bit harder for me to grasp ASPECT concepts 
because of the introspection and etc.).
However, I seem to have failed miserably at one point.
My understanding is that I can get
   A x = f- Pf

By:
- First I calculate my projection operator. It is the integral of the 
pressure shape function over ALL elements associated with a DOF. I 
calculate it by looping over the elements and integrating with a quadrature 
(see code below)
- Then I assemble my system and my right-hand side as usual
- Before I solve my linear system, I modify my Right-hand-side so that it 
does not include the constant mode. I litteraly calculate f= f- Pf
- I solve my linear system for this new RHS.

However, my GMRES stops very quickly after a certain number of newton 
iteration with the following : AztecOO::Iterate error code -4: GMRES 
Hessenberg ill-conditioned
Clearly I am doing something wrong / I understood something wrong regarding 
this, because I know the rest of my code works.
In my system there is no block matrices and everything is lumped in the 
same matrix. I havent re-ordered by components either.

Calculation of the projection:
template 
void
GLSNavierStokesSolver::initializePressureRHSCorrection()
{
  TimerOutput::Scope t(computing_timer, "pressure_RHS_projector");
  pressure_shape_function_integrals=0;
  QGauss   quadrature_formula(degreeQuadrature_);
  FEValues fe_values (fe,
   quadrature_formula,
   update_values |
   update_quadrature_points |
   update_JxW_values );

  const unsigned intdofs_per_cell = fe.dofs_per_cell;
  const unsigned intn_q_points= quadrature_formula.
size();
  const FEValuesExtractors::Scalar  pressure (dim);
  Vector   
 local_pressure_shape_function_integrals(dofs_per_cell);
  std::vector  local_dof_indices (dofs_per_cell);


  typename DoFHandler::active_cell_iterator
  cell = dof_handler.begin_active(),
  endc = dof_handler.end();
  for (; cell!=endc; ++cell)
{
  if (cell->is_locally_owned())
{
  fe_values.reinit(cell);
  local_pressure_shape_function_integrals=0;
  for (unsigned int i=0; iget_dof_indices (local_dof_indices);
  zero_constraints.distribute_local_to_global(
local_pressure_shape_function_integrals,
  local_dof_indices,
  
pressure_shape_function_integrals);
}
}
  pressure_shape_function_integrals.compress (VectorOperation::add);
}


Projection of the System RHS onto it's space without the constant:
template 
void GLSNavierStokesSolver::correctRHSClosedSystem()
{


  // calculate projection of system_rhs on pressure_shape_function_integrals
  for (unsigned int dof_i=0 ; dof_i < pressure_shape_function_integrals.size
() ; ++dof_i)
pressure_projection[dof_i] = pressure_shape_function_integrals[dof_i]*
system_rhs[dof_i];

  // Add it to RHS
  system_rhs.add(-1.,pressure_projection);
}


Thank you again for everything,
This is greatly appreciated.
Bruno


On Tuesday, 5 March 2019 16:10:09 UTC-5, Wolfgang Bangerth wrote:
>
>
> > A quick question. I think I understand what is done in 
> > 
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
>  
> > Weirdfully, I found the "pickaxe" version easier to understand (Love the 
> > comments by the way, this is awesome) 
>
> :-) I suspect that Timo gets credit for that one! 
>
>
> > A question I have is related to the _pressure_shape_function_integrals_ 
> member. 
> > My understanding is that this is the integration over the element cellI 
> of the 
> > shape function associated with pressure. 
>
> I believe it's actually the integral over all cells the shape function is 
> associated with. 
>
>
> > However, I have looked over the ASPECT documentation and I have not been 
> able 
> > to find where this is implemented. Is it a vector that is manually 
> filled or 
> > is there a helper function that can be used to automatically generate 
> it? 
>
> Yes, that required a bit of searching. This is filled here: 
>
> https://github.com/geodynamics/aspect/blob/master/source/simulator/assemblers/stokes.cc#L601
>  
>
> 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 

Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

Dear Wolfgang,
A quick question. I think I understand what is done in 
https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
Weirdfully, I found the "pickaxe" version easier to understand (Love the 
comments by the way, this is awesome)

A question I have is related to the *pressure_shape_function_integrals* 
member.
My understanding is that this is the integration over the element cellI of 
the shape function associated with pressure.
However, I have looked over the ASPECT documentation and I have not been 
able to find where this is implemented. Is it a vector that is manually 
filled or is there a helper function that can be used to automatically 
generate it?

Thanks for everything
Bruno

-- 
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] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

Dear Bruno and Wolfgang,

Thank you for your answers.
 I believe Wolfgang's answer is exactly what i had in mind (but said in 
clear words...). I will look at the Aspect code and try to port that to 
mine.
Thank you for the very detailed answer.
Best
Bruno


On Tuesday, 5 March 2019 10:08:29 UTC-5, Wolfgang Bangerth wrote:
>
>
> > One of the issue of my system is that the pressure is defined up to a 
> > constant. On coarse mesh this does not affect the GMRES solver. However, 
> on 
> > finer mesh, it seems that the GMRES Solver is greatly affected by this 
> > near-singularity of the matrix system. 
> > I have often read in the literature that for stabilized method, the best 
> way 
> > was to remove the "mode" associated to a constant pressure. I believe 
> this 
> > implies some sort of projection of the residual in a space without a 
> pressure 
> > constant? 
> There are two parts of this problem: 
>
> 1/ A deep theorem in linear algebra states that because the constant 
> pressures 
> are in the kernel of the matrix (i.e., Ay=0 for all vectors y that 
> correspond 
> to a constant pressure and zero velocity), that the *range* of the matrix 
> A 
> has dimension at most n-1. As a consequence, if you have a linear system 
>A x = f 
> then it is only possible to find a vector x with 
>|| A x - f ||  =  0 
> if f is in the range of the matrix A. That will not generally be the case, 
> due 
> to discretization and integration errors. If f is not in the range of A, 
> there 
> is no way for GMRES to reduce the residual below a certain threshold, no 
> matter how many iterations one runs. 
>
> The way to solve this is to solve 
>A x = f - Pf 
> instead where Pf is the projection onto the subspace not reachable by A. 
> This 
> is easy enough if you have a uniform mesh, but it's a bit complicated if 
> that's not the case. It's also complicated if one uses an enriched 
> pressure 
> space. Here is ASPECT's implementation of this step: 
>
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
>  
> You will be able to copy this and simplify it for your case. 
>
>
> 2/ There is now a null space and GMRES will find one of the infinitely 
> many 
> solutions of 
>A x = f 
> This may not be the solution you are looking for, so you probably want to 
> update the additive constant in the pressure after solution. How you want 
> to 
> do this depends on the application. Here again is the implementation in 
> ASPECT: 
>
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L753
>  
>
> I hope this helps! 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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)

[Bruno Blais]

Hello everyone,
I am currently solving a GLS stabilized form of the Navier-Stokes equation 
using DEALII.
The residual of the system looks similar to the regular Incompressible 
Navier-Stokes, except that a stiffness matrix that is dependent on the 
element size is added to the P-P block.
I have a great success / scaling when solving Manufactured Solutions and I 
have been able to reach relatively decent P1-P! meshes (say 2056x2056 or 
256x256x256) on a workstation computer.
I am currently using the GMRES Trilinos Wrapper with ILU preconditioning 
with a relatively high fill-in level (4).

One of the issue of my system is that the pressure is defined up to a 
constant. On coarse mesh this does not affect the GMRES solver. However, on 
finer mesh, it seems that the GMRES Solver is greatly affected by this 
near-singularity of the matrix system.
I have often read in the literature that for stabilized method, the best 
way was to remove the "mode" associated to a constant pressure. I believe 
this implies some sort of projection of the residual in a space without a 
pressure constant?

I know this is a very broad question, but how would I go and implement such 
a thing in my solver?
Are there any examples in DEALII where a Poisson problem is solved but with 
strictly Neumann boundaries? That would be a very similar problem that 
could guide me in the right direction.

Thank you very much for your help. The support of this forum is greatly 
appreciated. 

-- 
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] Shared memory parallelism of Trilinos

[Bruno Blais]

Dear Bruno,
Thanks, that definitely answers my questions.
I would have a final slightly related question. 
Would I see significant performance gain for a GMRES + ILU preconditioning 
setup by going from AztecOO to the Tpetra or Belos stack?
I guess that the best parallel performance increase I could get would be by 
using the AMG from ML, but I think this will require some experiments 
before I can set-up the parameters correctly. My system is not elliptic and 
can be quite ugly at times due to stabilization.

Thank you for everything.

On Monday, 4 March 2019 09:03:24 UTC-5, Bruno Turcksin wrote:
>
> Le lun. 4 mars 2019 à 08:44, Bruno Blais > 
> a écrit : 
> > I'm using the wrapper, so I guess by default that means it is using the 
> AztecOO stack of solvers? 
>
> Yes, that's right. You won't get any speedup using OpenMP with 
> AztecOO, you need to switch to the Tpetra stack and Belos to use 
> OpenMP (but we don't have wrappers for the all Tpetra stack) 
>
> >>   2) Why do you think that OpenMP would be faster than MPI? MPI is 
> usually faster than OpenMP unless you are very careful about your data 
> management. 
> > My original idea was that since in shared memory parallelism you could 
> precondition a larger chunk of the matrix as a whole, that the ILU 
> preconditioning would be more efficient in a shared-memory context than in 
> a distributed one. Thus you would need less GMRES iterations to solve > 
> your system. It seems I am wrong :) ? 
>
> Using larger blocks for ILU preconditioning will decrease the number 
> of GMRES iterations but you will spend more time in ILU, so it's hard 
> to say if it's worth it. 
>
> Best, 
>
> Bruno 
>

-- 
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] Shared memory parallelism of Trilinos

[Bruno Blais]

Bruno 

On Monday, 4 March 2019 08:28:18 UTC-5, Bruno Turcksin wrote:

> Bruno,
>
> On Monday, March 4, 2019 at 7:41:27 AM UTC-5, Bruno Blais wrote:
>>
>> 2. Furthermore, when you compile Trilinos with OpenMP and you try to 
>> compile the latest version of DEAL.II, you get a compilation error when 
>> ".hpp" from Kokkos are included. The error reads something like:
>> Kokkos was compiled with Openmp but the compiler did not pass an openmp 
>> flag.
>>
>> This can be easily fixed by manually adding -fopenmp to the CXX flags 
>> used by dealii. However, would it not be a better idea to add a 
>> DEALII_ENABLE_OpenMP Flag directly in the cmake to ensure that if you put 
>> that flag on, the -fopenmp flag is enabled?
>> Maybe I missed such an option. It just made me unsure if I was doing 
>> something supported or not.
>>
> deal.II does not use OpenMP for multithreading so if -fopenmp is missing 
> it's because Trilinos did not export it correctly. Unless you mean that you 
> include Kokkos in your own code? In that case you are responsible for the 
> flags if you use OpenMP. 
>

The issue would be then that my Trilinos installation did not export the 
flag correctly. I will check that out. It makes sense. I only use the 
solvers with the DEALII TrilinosWrappers, I have not tried to use the 
Trilinos solvers directly.
 

>
>> 3. When compiled with OpenMP, I got deceptively poor performances, but 
>> maybe this is because of the relatively small size of my application. I 
>> would have (naively maybe) expected that the time to solver a linear system 
>> with GMRES using 1 MPI with 4 OMP threads would have been lower than the 
>> time it takes with 4 MPI and on my application this was not the case. I was 
>> surprised because I was expecting my ILU Preconditioning to work better on 
>> lower amount of cores, but maybe this is related to fill-in or other issues?
>>
> Two things here:
>   1) Which package are using? The Epetra stack does not support OpenMP so 
> you can compile with OpenMP but it won't be used.
>

I'm using the wrapper, so I guess by default that means it is using the 
AztecOO stack of solvers?
 

>   2) Why do you think that OpenMP would be faster than MPI? MPI is usually 
> faster than OpenMP unless you are very careful about your data management.
>

My original idea was that since in shared memory parallelism you could 
precondition a larger chunk of the matrix as a whole, that the ILU 
preconditioning would be more efficient in a shared-memory context than in 
a distributed one. Thus you would need less GMRES iterations to solve your 
system. It seems I am wrong :) ?

 

>
> Best,
>
Thanks, this is very interesting / enlightening
 

>
> Bruno
>

-- 
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] Shared memory parallelism of Trilinos

[Bruno Blais]

Hello everyone,
I have a quick question for which I have not found documentation.
The suggested compilation options for Trilinos do not suggest the use of 
OpenMP and the flag is not enabled by default. 
DEALII by default also compiles using TBB for shared memory parallelism. 
Does this mean that the compilation of Trilinos with the suggested flag 
enables shared memory parallelism when using the TrilinosWrappers function, 
or is it necessary to compile Trilinos with openmp or threads enabled?
What are the official guidelines for this?

Thanks a lot
Bruno

-- 
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] Setting constant modes in Trillinos AMG preconditiner for (u,v,w,p) type of problems

[Bruno Blais]

Hello everyone,
I am currently solving a stabilized form of the Navier-Stokes equations. 
Due to the stabilization terms, I can easily lump everything in a single 
matrix.
Everything is done via the Trillinos Wrappers. 
I get relatively nice scaling up to a certain number of cells (say 1e6 or 
so), but then the resolution time greatly increases, especially in parallel.
My feeling is that this is due to the decreased efficiency of the ILU 
preconditioner that I use.
I would like to test the use of the AMG preconditioner of Trillinos.
I have managed to set-up everything accordling, but I have an issue with 
the constant mode in the additional data.

https://www.dealii.org/8.5.0/doxygen/deal.II/structTrilinosWrappers_1_1PreconditionAMG_1_1AdditionalData.html#af4c9b8fcda773646bafc45a09b9c800f

Which reads:

*Specifies the constant modes (near null space) of the matrix. This 
parameter tells AMG whether we work on a scalar equation (where the near 
null space only consists of ones, and default value is OK) or on a 
vector-valued equation. For vector-valued equation problem with 
n_component, the provided constant_modes should fulfill the following 
requirements: *

   - * n_component.size() == n_component* 
   - * n_component[*].size() == n_dof_local or n_component[*].size() == 
   n_dof_global *
   - * n_component[ic][id] == "idth DoF is corresponding to component ic* 

My understanding is that since my problem has 4 DOF per node (if I am P1), 
I need to use specify the correct constant_modes. However, I have not 
really been able to figure out how. I know this may sound like a very 
simple question, but could someone guide me in how I should specify by 
constant modes if I have a nDim+1 component problem? I don't think I really 
understand this part of the documentation.

Thank you for your understanding and help.

-- 
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] Choosing the appropriate parameters for TrilinosWrappers::PreconditionILU for GLS stabilized Navier-Stokes

[Bruno Blais]

Hello everyone,
I am currently solving a stabilized version of the Navier-Stokes equation 
using dealii and everything is working well when I use a direct solver.
I would like (obviously?) to extend it to use a FGMRES iterative solver to 
be able to tackle very large systems. I have ported those aspects to use 
Trilinos.

The weak form  is (for now we can neglect the temperature related term) :


Concretely, the additional term due to the stabilization affect the 
velocity block in a SUPG manner and the pressure block by adding an 
element-dependent stiffness matrix.
I have read that ILU and ILUT type of preconditioners are the best-suited  
for this type of systems, but I find myself unable to orient my choice of 
parameters.
When we look at the constructor for the AdditionalData of ILU and ILUT 
preconditioner there are the following parameters:
ILU :
  AdditionalData(const unsigned int ilu_fill = 0,
 const double   ilu_atol = 0.,
 const double   ilu_rtol = 1.,
 const unsigned int overlap  = 0);
ILUT:
  AdditionalData(const double   ilut_drop = 0.,
 const unsigned int ilut_fill = 0,
 const double   ilut_atol = 0.,
 const double   ilut_rtol = 1.,
 const unsigned int overlap   = 0);

The documentation suggest an atol between 1e-2 and 1e-5 and a rtol of the 
order of 1.01. However, what should my ilu_fill and ilut_drop be? I 
understand that increasing the fill will increase the memory consumption 
and that I should try to reach a compromise, but right now I find that it's 
mostly the speed of my iterative solve that is the problem. Are there any 
guidelines to orient my choice of parameters? Am I going with the right 
type of pre-conditioner or am I in the wrong here?
Thank you very much for your help.

-- 
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] Change in the behavior of Boundary Id interpretation from GMSH 2 to GMSH 4?

[Bruno Blais]

Dear Daniel,
Sorry I dropped the ball on this one.
I checked out the patch and tested it on two of my full cases and 
everything works well.
No differences between a GMSH ASCII 2 and GMSH ASCII 4, whereas before I 
had the previous issue.

>From my point of view this fixes everything.

Thank you very much!

On Tuesday, 29 January 2019 20:05:36 UTC-5, Daniel Arndt wrote:
>
> Bruno,
>
> Can you try if the patch in https://github.com/dealii/dealii/pull/7637 
> fixes your problem?
>
> 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.

Re: [deal.II] Re: the viscous term in SUPG stabilization terms

[Bruno Blais]

I might be late to the party,
But I believe what Wolfgang is trying to say is that Q1-Q1 elements are 
bi-linear. Consequently, except if you mesh is perfectly alligned with the 
x and y axis, you will have a non-zero second derivative of the shape 
function.
Your shape function are of the form : a + b*x + c*y + d*x*y in your 
reference element. Thus if your Q1-Q1 element is slightly tilted you will 
get a non-zero second derivative.


On Wednesday, 17 January 2018 16:13:42 UTC-5, Feimi Yu wrote:
>
> Thanks a lot, Daniel and Wolfgang!
>
> What I understood was that since the SUPG and PSPG terms are integrated 
> over the cell interiors so the laplacians can be ignored while using Q1/Q1 
> element without any singularity problems on the cell edges.
>
> Thanks!
> Feimi
>
> On Wednesday, January 17, 2018 at 2:36:39 PM UTC-5, Wolfgang Bangerth 
> wrote:
>>
>> On 01/17/2018 11:59 AM, Daniel Arndt wrote: 
>> > These terms are simply zero if you consider Q1 elements for the 
>> velocity. I 
>>
>> That's true only if you're on affine meshes. In general, the Laplacian 
>> of a Q1 function on an arbitrary mesh is not zero. 
>>
>> That doesn't mean, though, that it may not be valid to just ignore the 
>> issue. 
>>
>> But I think what you want to say is that the term is supposed to be 
>> interpreted as a sum over cell interiors, not as an integral over the 
>> entire domain. In the latter case, you'd have to worry about the fact 
>> that Delta uh is singular on edges, whereas in the former you ignore 
>> these contributions. 
>>
>> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Change in the behavior of Boundary Id interpretation from GMSH 2 to GMSH 4?

[Bruno Blais]

Daniel,
Thank you very much. I real the issue and saw that you guys found a smaller 
mesh to break the example.
I'll try to look further into it in the following weeks and see if I can 
come up with a nice solution.
Best
Bruno


On Tuesday, December 4, 2018 at 9:04:13 AM UTC-5, Daniel Arndt wrote:
>
> Bruno,
>
> I created a corresponding issue in the GitHub repository at 
> https://github.com/dealii/dealii/issues/7501.
>
> 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.

Re: [deal.II] Change in the behavior of Boundary Id interpretation from GMSH 2 to GMSH 4?

[Bruno Blais]

Dear Wolfgang,
Thanks for looking into it.
I can't guarantee a timeline, but I'll try to take a look at it this 
weekend / early next week.
I'll use  your code to debug, it is indeed significantly smaller!. I'll try 
to see if I can come up with a very small mesh that breaks it down (like 
10cells or so) in a similar fashion and try to work through it afterward.

Thanks for the help, i'll keep you posted :)!
Best,
Bruno


On Thursday, 29 November 2018 00:11:11 UTC-5, Wolfgang Bangerth wrote:
>
>
> Bruno, 
>
> > I have been having some issues when trying to run a "toy" code (a simple 
> > Steady Navier-Stokes solver) on a case I used to run. 
> > The boundary group Id 0, which is manually set in my GMSH script 
> > (vonKarman.geo), is detected perfectly as both a BC and a manifold when 
> I 
> > export to GMSH 2 and I get the expected results. 
> > However, if I export the same GMSH file to V4 the same code does not run 
> > anymore. I guess there is a change in behavior I did not notice in how 
> dealii 
> > interprets the new format? What am I doing wrong? 
> > The error I get is (when running in debug mode): 
> > | 
> >  
> > Anerror occurred inline <10387>of file 
> > infunction 
> > 
> voiddealii::Triangulation::set_all_manifold_ids_on_boundary(dealii::types::boundary_id,dealii::types::manifold_id)[withintdim
>  
>
> > =2;intspacedim =2;dealii::types::boundary_id 
> > =unsignedint;dealii::types::manifold_id =unsignedint] 
> > Theviolated condition was: 
> >  boundary_found 
> > Additionalinformation: 
> > Thegiven boundary_id 0isnotdefinedinthisTriangulation! 
> > 
> > Stacktrace: 
> > --- 
> > #0  /home/blaisb/codes/dealii/build/lib/libdeal_II.g.so.9.1.0-pre: 
> > dealii::Triangulation<2, 2>::set_all_manifold_ids_on_boundary(unsigned 
> int, 
> > unsigned int) 
> > #1  ./schursteadynavierstokes: 
> SchurNavierStokesSolver<2>::runVonKarman() 
> > #2  ./schursteadynavierstokes: main 
> >  
> > 
> > | 
>
> I put this piece of code into your source file right after reading in the 
> mesh: 
> | 
>  for (const auto  : triangulation.active_cell_iterators()) 
>for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; 
> ++f) 
>  if (cell->face(f)->at_boundary()) 
>std::cout << "cell=" << cell << ", face=" << f 
>  << ", boundary_id=" << cell->face(f)->boundary_id() 
>  << std::endl; 
> | 
> The output shows that all boundary faces have boundary ids between 1 and 
> 8. 
>
> Of course you already know that from the error message. The question is 
> how 
> these get there. I have to say that I don't understand the gmsh file 
> format 
> version 4 -- it's quite complicated with the entities. But do you think 
> you 
> could try to debug this problem by stepping into the GridIn and seeing 
> where 
> and why the SubCellData objects are created that assign nonzero 
> boundary_ids 
> to the boundary faces? The gmsh format 4 reader is new (just a few weeks) 
> and 
> it would not be a surprise if it still had bugs in it. We would definitely 
> appreciate any help with this! 
>
> A good step would also be to create a smaller mesh that illustrates the 
> issue. 
> The mesh you have has 400+ cells and 17 entities. It is difficult to read 
> through the mesh file and understand what each part does. 
>
> For all further work, the following, much-reduced program also shows the 
> issue: 
> | 
> #include  
> #include  
>
>
> int main () 
> { 
>using namespace dealii; 
>
>const unsigned int dim = 2; 
>
>Triangulation<2> triangulation; 
>
>  GridIn grid_in; 
>  grid_in.attach_triangulation (triangulation); 
>  std::ifstream input_file("vonKarman.msh"); 
>  grid_in.read_msh(input_file); 
>
>  for (const auto  : triangulation.active_cell_iterators()) 
>for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; 
> ++f) 
>  if (cell->face(f)->at_boundary()) 
>std::cout << "cell=" << cell << ", face=" << f 
>  << ", boundary_id=" << cell->face(f)->boundary_id() 
>  << std::endl; 
>
>
>  triangulation.set_all_manifold_ids_on_boundary(0,0); 
> } 
> | 
>
> Best 
>   W. 
>
>
> > However, BC ID 0 as a physical group is defined explicitely in the 
> script. I 
> > have joined the .h and .cc of my code ( I am sorry it is ugly as hell, 
> it was 
> > just made for a quick test). 
> > 
> > The GMSH file is as follow (also in the tar archive): 
> > | 
> > 
> > y0=0; 
> > y1=1; 
> > 
> > x0 =0.; 
> > x1 =10; 
> > 
> > xc=1; 
> > yc=0.5; 
> > r=0.025; 
> > 
> > lo =2.0e-1; 
> > lc =1.0e-1; 
> > 
> > Point(0)={x0,y0,0,lo}; 
> > Point(1)={x0,y1,0,lo}; 
> > Point(2)={x1,y0,0,lo}; 
> > Point(3)={x1,y1,0,lo}; 
> > 
> > Point(4)={xc,yc,0,lc}; 
> > Point(5)={xc-r,yc,0,lc}; 
> > Point(6)={xc,yc+r,0,lc}; 
> > Point(7)={xc+r,yc,0,lc}; 
> > 

[deal.II] Change in the behavior of Boundary Id interpretation from GMSH 2 to GMSH 4?

[Bruno Blais]

Hello,
I have been having some issues when trying to run a "toy" code (a simple 
Steady Navier-Stokes solver) on a case I used to run.
The boundary group Id 0, which is manually set in my GMSH script 
(vonKarman.geo), is detected perfectly as both a BC and a manifold when I 
export to GMSH 2 and I get the expected results.
However, if I export the same GMSH file to V4 the same code does not run 
anymore. I guess there is a change in behavior I did not notice in how 
dealii interprets the new format? What am I doing wrong?
The error I get is (when running in debug mode):

An error occurred in line <10387> of file  in function
void dealii::Triangulation::
set_all_manifold_ids_on_boundary(dealii::types::boundary_id, dealii::types::
manifold_id) [with int dim = 2; int spacedim = 2; dealii::types::boundary_id 
= unsigned int; dealii::types::manifold_id = unsigned int]
The violated condition was: 
boundary_found
Additional information: 
The given boundary_id 0 is not defined in this Triangulation!

Stacktrace:
---
#0  /home/blaisb/codes/dealii/build/lib/libdeal_II.g.so.9.1.0-pre: 
dealii::Triangulation<2, 2>::set_all_manifold_ids_on_boundary(unsigned int, 
unsigned int)
#1  ./schursteadynavierstokes: SchurNavierStokesSolver<2>::runVonKarman()
#2  ./schursteadynavierstokes: main



However, BC ID 0 as a physical group is defined explicitely in the script. 
I have joined the .h and .cc of my code ( I am sorry it is ugly as hell, it 
was just made for a quick test).

The GMSH file is as follow (also in the tar archive):

y0=0;
y1=1;

x0 =0.;
x1 =10;

xc=1;
yc=0.5;
r=0.025;

lo = 2.0e-1;
lc = 1.0e-1;

Point(0) = {x0, y0, 0, lo};
Point(1) = {x0, y1, 0, lo};
Point(2) = {x1, y0, 0, lo};
Point(3) = {x1, y1, 0, lo};

Point(4) = {xc, yc, 0, lc};
Point(5) = {xc-r, yc, 0, lc};
Point(6) = {xc, yc+r, 0, lc};
Point(7) = {xc+r, yc, 0, lc};
Point(8) = {xc, yc-r, 0, lc};


// Contour Box
Line(1) = {0,1};
Line(2) = {1,3};
Line(3) = {3,2};
Line(4) = {2,0};

// Inner Circle
Circle(5)={5,4,6};
Circle(6)={6,4,7};
Circle(7)={7,4,8};
Circle(8)={8,4,5};

Line Loop(1) = {1,2,3,4};
Line Loop(2) = {5,6,7,8};

// Creates the physical entities
Physical Line(0)={5,6,7,8};
Physical Line(1)={1};
Physical Line(2)={2,4};
Physical Line(3)={3};

Plane Surface(1) = {1,2};
Physical Surface(2) = {1};
Recombine Surface{1};



The DEALII lines regarding the manifold are :
grid_in.read_msh(input_file);
Point circleCenter(1.0,0.5);
static const SphericalManifold boundary(circleCenter);
triangulation.set_all_manifold_ids_on_boundary(0,0);
triangulation.set_manifold (0, boundary);

And the lines linked to the BC are :
emplate 
void SchurNavierStokesSolver::setup_dofs ()
{
system_matrix.clear();
pressure_mass_matrix.clear();

dof_handler.distribute_dofs(fe);

std::vector block_component(dim+1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise (dof_handler, block_component);
dofs_per_block.resize (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, 
block_component);
unsigned int dof_u = dofs_per_block[0];
unsigned int dof_p = dofs_per_block[1];

FEValuesExtractors::Vector velocities(0);
{
  nonzero_constraints.clear();
  DoFTools::make_hanging_node_constraints(dof_handler, 
nonzero_constraints);
  VectorTools::interpolate_boundary_values(dof_handler,
   0,
   ZeroFunction(dim+1),
   nonzero_constraints,
   fe.component_mask(velocities
));

  if (simulationCase_==TaylorCouette)
  {
  VectorTools::interpolate_boundary_values(dof_handler,
   1,
   RotatingWall(),
   nonzero_constraints,
   fe.component_mask(
velocities));
  }
  if (simulationCase_==BackwardStep)
  {
  VectorTools::interpolate_boundary_values(dof_handler,
   1,
   PoiseuilleInlet(),
   nonzero_constraints,
   fe.component_mask(
velocities));
  }



  if (simulationCase_==vonKarman)
  {
  std::set no_normal_flux_boundaries;
  no_normal_flux_boundaries.insert (2);
  VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
  
 no_normal_flux_boundaries,
  
 nonzero_constraints
  

Re: [deal.II] Issue with convergence of iterative linear solver for system matrix in modified step-57 with no-normal flux constraints

[Bruno Blais]

I have found my problem. A boundary condition in my mesh was ill-defined 
and this lead to my error.
Sorry for the confusion.
Thanks


On Monday, 19 February 2018 13:41:35 UTC-5, Bruno Blais wrote:
>
> Hello,
> Sorry I feel I have not explained myself correctly. Here is a drawing of 
> the case:
>
>
> <https://lh3.googleusercontent.com/-4py-EISXcqc/WosZebbKS-I/C5w/-qHRMSkwW6IinnVGocwlonnFHrF3rXXbQCLcBGAs/s1600/text4174-2.png>
>
> With Ux may be a profile or a constant.
> Initially I had set-up that case with Ux as a Poiseuille flow profile and 
> the Top and Bottom as no slip boundary conditions. In that case I got the 
> solution I was expecting. Convergence was good and the solution was what I 
> had in mind.
> Then what I wish to do is set Ux as a constant and apply slip at the top 
> and bottom boundary conditions. When I do this is when I am unable to reach 
> convergence in my iterative solution, yet I do not understand why because 
> the problem should be well posed.
> I don't think that leaving my right boundary as a free traction boundary 
> condition is problematic since it is by all mean an outlet boundary. This 
> worked well for instance in a backward facing step case.
> Thank you for your time, sorry for the confusion.
> Best regards
>
>
>
> On Monday, 19 February 2018 10:50:48 UTC-5, Timo Heister wrote:
>>
>> > Does the order in which I apply the nonzero and zero constraints 
>> matter? 
>>
>> These are two independent objects, so no. 
>>
>> > Currently I apply the inlet and then the no-slip in the 
>> nonzero_constraints, 
>> > thus the bottom and top wall appear after the inlet. Afterward the 
>> cylinder 
>> > is put in the zero constraints. 
>>
>> I don't understand. You need boundary conditions in 
>> nonzero_constraints and zero_constraints for all boundaries. The 
>> reason for these two sets is that one is used for the initial solve, 
>> while the second one is used for the Newton updates. 
>>
>>
>> -- 
>> Timo Heister 
>> http://www.math.clemson.edu/~heister/ 
>>
>

-- 
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] Issue with convergence of iterative linear solver for system matrix in modified step-57 with no-normal flux constraints

[Bruno Blais]

Hello,
Sorry I feel I have not explained myself correctly. Here is a drawing of 
the case:



With Ux may be a profile or a constant.
Initially I had set-up that case with Ux as a Poiseuille flow profile and 
the Top and Bottom as no slip boundary conditions. In that case I got the 
solution I was expecting. Convergence was good and the solution was what I 
had in mind.
Then what I wish to do is set Ux as a constant and apply slip at the top 
and bottom boundary conditions. When I do this is when I am unable to reach 
convergence in my iterative solution, yet I do not understand why because 
the problem should be well posed.
I don't think that leaving my right boundary as a free traction boundary 
condition is problematic since it is by all mean an outlet boundary. This 
worked well for instance in a backward facing step case.
Thank you for your time, sorry for the confusion.
Best regards



On Monday, 19 February 2018 10:50:48 UTC-5, Timo Heister wrote:
>
> > Does the order in which I apply the nonzero and zero constraints matter? 
>
> These are two independent objects, so no. 
>
> > Currently I apply the inlet and then the no-slip in the 
> nonzero_constraints, 
> > thus the bottom and top wall appear after the inlet. Afterward the 
> cylinder 
> > is put in the zero constraints. 
>
> I don't understand. You need boundary conditions in 
> nonzero_constraints and zero_constraints for all boundaries. The 
> reason for these two sets is that one is used for the initial solve, 
> while the second one is used for the Newton updates. 
>
>
> -- 
> Timo Heister 
> http://www.math.clemson.edu/~heister/ 
>

-- 
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] Issue with convergence of iterative linear solver for system matrix in modified step-57 with no-normal flux constraints

[Bruno Blais]

Hello,
My gmsh mesh produces 3 physical groups which are independent. Say my 
simulation is the flow around a cylinder. The cylinder is a physical group 
(0), the inlet (say left wall) is a class inheriting from the Function 
class which sets the X component to 1 and the Y component to 0 and uses the 
velocity mask.
Finally the physical group(2) is supposed to be the top and the bottom 
boundary which would be slip. The last part (the right wall) I wish to 
leave as a free boundary condition so it is not put in any physical group.
Everything works if I put the physical group(2) as a dirichlet boundary 
condition (and if I adapt the velocity profile of the inlet so that it 
matches the 0 on the top and bottom wall using for instance a parabolic 
inlet profile).

Does the order in which I apply the nonzero and zero constraints matter?
Currently I apply the inlet and then the no-slip in the 
nonzero_constraints, thus the bottom and top wall appear after the inlet. 
Afterward the cylinder is put in the zero constraints.
Thanks you very much
Best regards
Bruno


On Sunday, 18 February 2018 11:58:40 UTC-5, Timo Heister wrote:
>
> Are you setting the conditions correctly for the zero_constraints and 
> the nonzero_constraints and are you sure you are not applying other 
> boundary conditions on that wall (for example with 
> interpolate_boundary_values())? 
>
> On Sat, Feb 17, 2018 at 5:46 PM, Bruno Blais <blais...@gmail.com 
> > wrote: 
> > Hello everyone, 
> > 
> > I believe this may sound like a relatively dumb question, but I thank 
> you 
> > for your time. 
> > I am using a (slightly) modified version of step-57 to solve certain 
> steady 
> > state Navier-Stokes problem. 
> > I have had relatively good success and showed that i could recover the 
> > appropriate order of convergence on manufactured solutions and obtained 
> good 
> > solutions on problems like a backward facing step, as such I am 
> confident 
> > that the residual / jacobian matrix / linear solver / Schur complement 
> > aspect is ok. Anyway, I did not modify anything from Step-57 when it 
> comes 
> > to the Schur complement / solution of the linear system. 
> > 
> > My core issue arises when I try to replace one of my homogenous or 
> > non-homogenous dirichlet boundary condition with a 
> > no_normal_flux_constraints to impose slip instead of no-slip on a 
> boundary. 
> > Simply, I can say that I implement it by adding an additional constraint 
> in 
> > the setup_dofs member function of step-57 such as: 
> > 
> > std::set no_normal_flux_boundaries; 
> > 
> > no_normal_flux_boundaries.insert (2); /* here 2 is a Physical Line in my 
> > gmsh mesh */ 
> > 
> > VectorTools::compute_no_normal_flux_constraints (dof_handler, 0, 
> > 
> > 
>  no_normal_flux_boundaries, 
> > 
> >  nonzero_constraints 
> > 
> >  ); 
> > 
> > 
> > Now everything compiles , but the iterative solver (GMRES in this case) 
> for 
> > the system matrix does not converge anymore with this set of boundary 
> > conditions. My general problem is well-posed and if I replace these 
> boundary 
> > with regular Dirichlet I get an expected solution. 
> > Clearly, I am doing something wrong, but I must admit my lack of 
> knowledge 
> > on the issue. Could it be related to the way for system matrix and the 
> Schur 
> > complement is formed which is rendered not-ok in this case for no-flux 
> > boundaries? 
> > I know the example is made with Dirichlet (homogenous or not) boundary 
> in 
> > mind, so clearly I am missing something when it comes to changing one of 
> > these boundaries to a slip. 
> > 
> > I thank you greatly for your time, 
> > Best regards, 
> > Bruno 
> > 
> > 
> > -- 
> > The deal.II project is located at 
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_=DwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=ulWuM0dNK_XHF-DL_x2EvT8bMSBuSrV4XKcqEfBiLm8=-4wmXq2eve0gPycm8fOaR7d-hpgfjIE1aPscZQsUILk=
>  
> > For mailing list/forum options, see 
> > 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_forum_dealii-3Fhl-3Den=DwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=ulWuM0dNK_XHF-DL_x2EvT8bMSBuSrV4XKcqEfBiLm8=shEU_TxDkAEzy7XeIgIrmbDaeclYOji85cAxehiBlJE=
>  
> > --- 
> > You received this message because you are subscribed to the Google 
>