[deal.II] Re: How to implement GridGenerator::channel_with_cylinder?

2019-01-15 Thread nboychenko1
Does anybody have an example of how to fill out the hole and assign to it 
different material ID?

Thank you,
Nick.

On Wednesday, December 19, 2018 at 12:07:16 PM UTC-5, Daniel Arndt wrote:
>
> Chucui,
>
> the mesh you are showing degenerates close to the midpoint. This isvery 
> likely not what you wan't.
> For the inner part something like shown in 
> https://www.dealii.org/developer/doxygen/deal.II/classTransfiniteInterpolationManifold.html
> would be a good idea.
> Then for the middle part GridGenerator::hyper_cube_with_cylindrical_hole (
> https://www.dealii.org/9.0.0/doxygen/deal.II/namespaceGridGenerator.html#add14cab546d033c1eaacc9234c64ebcd
> )
> should work, you just have to make sure that the vertices align when 
> merging the two meshes via 
> GridGenerator::merge_triangulations (
> https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#a7cd88e7eacd46697dee80ad2b8438d54
> ).
> Make sure too choose copy_manifold_ids if you are using a developer 
> version or set the manifold ids after creating and merging the whole mesh.
> Finally, for the outer part you would just merge some 
> GridGenerator::subdivided_hyper_rectangle (
> https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971
> ).
> Again make sure that the verices align and make sure the manifold ids are 
> correct after merging all the meshes.
>
> 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: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

2019-01-15 Thread Wolfgang Bangerth
On 1/15/19 10:02 AM, ntets...@berkeley.edu wrote:
> thank you for replying. I have looked thoroughly step 45 but it's not of much 
> help for my case. I actually want to impose periodic boundary conditions for 
> both velocities and pressure (I need to have a representative periodic cell).

Is this really what you want to do? Can you elaborate on your setup?

I'm asking this because generally, mathematically speaking, you can only 
impose the *stress* at the boundary, not the pressure alone. Second, if you 
think for example about an infinite pipe and you cut a piece out, then at the 
left and right end of that pipe segment, you can enforce that the velocity is 
periodic -- but the pressure will *not* be periodic, as there will be a 
pressure drop along the pipe which counters the friction forces.

Best
  W.

-- 

Wolfgang Bangerth  email: bange...@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] Re: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

2019-01-15 Thread ntetsika
Hi Daniel,

thank you for replying. I have looked thoroughly step 45 but it's not of 
much help for my case. I actually want to impose periodic boundary 
conditions for both velocities and pressure (I need to have a 
representative periodic cell). I used the following lines to set pressure 
to zero at the first node, but still gives me the same problem when I try 
to invert [B^T (diag M)^{-1} B] with cg solver for the purpose of Block 
Preconditioning - it says [B^T (diag M)^{-1} B] not invertible.

// constraint first pressure dof to be 0

std::vector boundary_dofs (dof_handler.n_dofs(), false);

DoFTools::extract_boundary_dofs (dof_handler, fe.component_mask(pressure), 
boundary_dofs);

const int first_boundary_dof = std::distance (boundary_dofs.begin(), 
std::find (boundary_dofs.begin(), boundary_dofs.end(), true));

constraints.add_line(first_boundary_dof);


I was thinking maybe the time-dependent stokes problem I'm trying to solve 
with fully periodic boundary conditions (velocities+pressure) is 
over-constrained and maybe that's why [B^T (diag M)^{-1} B] is not 
invertible. What do you think? Any ideas?


Thanks again for the help!


Best,

Magda

On Tuesday, January 15, 2019 at 11:05:38 AM UTC-5, Daniel Arndt wrote:
>
> Magda,
>
> [...]
>> DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, 
>> constraints);
>> DoFTools::make_periodicity_constraints(dof_handler, 4, 1, 0, 
>> constraints);
>>
>>
> These two lines imply that you want to impose periodic boundary conditions 
> for all the components. Assuming the dof_handler variable is of type 
> FESystem and describes both velocity and pressure (as is the case in the 
> code-gallery example),
> this in particular holds true for the pressure. Is that what you are 
> intending?
> If so, you probably need another constraint to make the pressure unique.
> Otherwise, you might want to have a look at step-45 
>  where a 
> Stokes problem is solved with periodic boundary conditions on part of the 
> boundary.
> There only periodicity constraints for the velocity are stored in the 
> ConstraintMatrix object.
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Re: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

2019-01-15 Thread Daniel Arndt
Magda,

[...]
> DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, constraints);
> DoFTools::make_periodicity_constraints(dof_handler, 4, 1, 0, constraints);
>
>
These two lines imply that you want to impose periodic boundary conditions 
for all the components. Assuming the dof_handler variable is of type 
FESystem and describes both velocity and pressure (as is the case in the 
code-gallery example),
this in particular holds true for the pressure. Is that what you are 
intending?
If so, you probably need another constraint to make the pressure unique.
Otherwise, you might want to have a look at step-45 
 where a Stokes 
problem is solved with periodic boundary conditions on part of the boundary.
There only periodicity constraints for the velocity are stored in the 
ConstraintMatrix object.

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] Several questions about mesh generation and distributed triangulations

2019-01-15 Thread Pascal Kraft
Dear Wolfgang,

thank you for your reply! I only noticed it now a while later since I had 
thought the topic was dead back when I posted it. Thank you for your time 
and effort!
Remarks to your response below:

Am Dienstag, 30. Oktober 2018 17:26:21 UTC+1 schrieb Wolfgang Bangerth:
>
>
> Pascal, 
> I don't think anyone responded to your email here: 
>
> > I will try to be as short as possible 
>
> That was only moderately successful ;-)) 
>
>
> > - if more details are required 
> > feel free to ask. Also I offer to submit all mesh generation code I 
> > create in the future, since others might have similar needs at some 
> point. 
>
> We would of course love to include this! 
>
I will, once I am satisfied with my solution, offer some examples and make 
them available. 

>
>
> > I work on a 3d mesh with purely axis-parallel edges. The mesh is a 
> > 2d-mesh (say in x and y direction) extruded to 3d (z-direction). Due to 
> > the scheme I use it is required, that the distributed version of this 
> > mesh be distributed as intervals along the z-axis ( process 0 has all 
> > cells with cell->center()[2] in [0,1], process 1 has all cells with 
> > cell->center()[2] in (1,2] and so on.) 
> > What I did originally was simply generating a mesh consisting of 
> > n_processes cells, let that mesh be auto partitioned, then turning 
> > partitioning off 
>
> Out of curiosity, how do you do this? 
>
I start with p::d::tria with 
parallel::distributed::Triangulation<3>::Settings::no_automatic_repartitioning  
and hand it to GridGenerator::subdivided_parallelepiped<3, 3>. In this call 
I use the version that also takes a vector of subdivision per by dimension 
and hand it a vector containtin [1,1,n_processes]. I then manually call 
tria->repartition() which gives me an ordered triangulation where the cells 
are partitioned according to their z-coordintates. I can then perform 
refinement and the basic structure remains the same because automatic 
repartitioning is turned off. The only real downside to this (apart from 
feeling really weird and hacky) is that I am now stuck with a p::d::tria 
which doesn't allow anisotropic refinement which would be nice for me and I 
cannot use another starting mesh then the mesh with one cell per processor 
because then the partitioning done by p4est is no longer easy to predict.

>
>
> > So my questions would be: 
> > 1. given a 2D-mesh, what is the easiest way to get a distributed 
> > 3d-extrusion with the partitioning described above?  (auto-partitioning 
> > generates balls in this mesh, not cuboids with z-orthogonal 
> > process-to-process interfaces) One suggestion for a function here would 
> > be a function to transform a shared to a distributed mesh because on a 
> > shared mesh I could set the subdomain Ids and then just call that 
> > function when I'm done 
>
> You can't do this for a parallel::distributed::Triangulation. That's 
> because the partitioning is done in p4est and p4est orders cells in the 
> depth-first tree ordering along a (space filling) curve and then just 
> subdivides it by however many processors you have. 
>
> The only control you have with p4est is how much weight each cell has in 
> this partition, but you can't say that cutting between cells A and B 
> should be considered substantially worse than cutting between cells C 
> and D -- which is what you are saying: Cutting cells into partitions is 
> totally ok if the partition boundary runs between copies of the base 2d 
> mesh, but should be prohibited if the cut is between cells within a copy. 
>
> Other partitioning algorithms allow this sort of thing. The partition 
> graphs (where each node is a cell, and an edge the connection between 
> cells). Their goal is to find partitions of roughly equal weight (where 
> the weight of a partition is the sum of its node weights) while trying 
> to minimize the cost (where the cost is the sum of the edge weights over 
> all cut edges). If you were to assign very small weights to all edges 
> between cells of different copies, and large weights to edges within a 
> copy, then you would get what you want. 
>
> It should not be terribly difficult to do this with 
> parallel::shared::Triangulation since there we use regular graph 
> partitioning algorithms. But I don't see how it is possible for 
> parallel::distributed::Triangulation. 
>
> OK, that is what I feared. p::s::T would work but I will keep looking for 
a solution with p::d::T since this application should scale to huge meshes 
eventually. 

>
> > 2. say I have a layer of faces in that mesh (in the interior) and I need 
> > nedelec elements and nodal elements on these faces (codimension 1) to 
> > evaluate their shape function and gradients in quadrature points of a 
> > 2D-quadrature on the faces. What is the best way to do this? (if it was 
> > a boundary I could give a boundary id and call 
> > GridGenerator::extract_boundary_mesh but it's not always on the 
> boundary. 
>
> I don't think it 

Re: [deal.II] Non-linear reaction diffusion PDE

2019-01-15 Thread Wolfgang Bangerth
On 1/14/19 5:56 PM, D. Sarkar wrote:
> 
> Thanks for the suggestions. Apologies for including the governing equation.
> The equation is shown below with a non-linear reaction term:
> 
> 
> 
> I had looked at the example you suggested but will have to go over your video 
> lectures.
> Using some assumptions, I was able to linearize the PDE and followed the 
> transient heat equation example (step-26) by absorbing the linearized 
> reaction 
> term with the spatiotemporal source term.

Yes, that's the way to do it -- treat the diffusion term with an implicit time 
discretization and make the reaction term explicit. You will have to choose a 
small time step this way, but the advantage is that the equation then 
essentially looks like the heat equation of step-26. If you want to choose 
larger time steps, you can either be semi-implicit, for example using

   a u^n / (b + u^{n-1})

for the reaction term, or make it fully implicit

   a u^n / (b + u^n)

in which case you end up with a nonlinear system to solve in each time step -- 
which you would then do using something like step-15.

Best
  W.

-- 

Wolfgang Bangerth  email: bange...@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] Problem with Block Preconditioner and Fully Periodic Boundary Conditions

2019-01-15 Thread ntetsika


Hi,


any chance I could get your help with that? I'm trying to use the code of 
@Jie-Cheng  for Block Preconditioner 
(https://github.com/dealii/code-gallery/tree/893bdb5b0a770758fce98859ab4d137f95e62a50/time_dependent_navier_stokes)
 
for a case of fully periodic boundary conditions, but then it gives me the 
error that [B^T (diag M)^{-1} B]^{-1} not invertible. Any ideas why this 
could happen or how could I fix it? I have increased the number of 
iterations for the CG solver and I also tried to use gmres, but none of 
them ways worked.


The way I impose the Periodic Boundary conditions is the following. I have 
a unit square domain with boundary ids 1,2,3 and 4 for the four sides.


std::vector > periodicity_vector;
std::vector > periodicity_vector1;

FullMatrix rotation_matrix(dim);
rotation_matrix[0][0]=1.;
rotation_matrix[1][1]=1.;

GridTools::collect_periodic_faces(triangulation, 2, 3, 1,
  periodicity_vector, Tensor<1, dim>(),
  rotation_matrix);
GridTools::collect_periodic_faces(triangulation, 4, 1, 0,
  periodicity_vector1, Tensor<1, dim>(),
  rotation_matrix);

triangulation.add_periodicity(periodicity_vector);
triangulation.add_periodicity(periodicity_vector1);

DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, constraints);
DoFTools::make_periodicity_constraints(dof_handler, 4, 1, 0, constraints);


Thank you for your help in advance, I've been stuck with that a long time 
now and I really don't understand what's the problem.


Best,
Magda

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


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

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

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

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


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2009 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

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


#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 

#include 
#include 

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

#include 
#include 

#include 
#include 
#include 

#include 
#include 


#include 
#include 


namespace Step37
{
  using namespace dealii;


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



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

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

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



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



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



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

LaplaceOperator();