[deal.II] Lecture 31.65(Picard iteration)

[hanks0227]

Hi all,

I have a question in the lecture 31.65 that is about Picard iteration.

In the lecture, the professor doesn't include   system_matrix.reinit 
(sparsity_pattern),  system_rhs.reinit (dof_handler.n_dofs()), and 
solution.reinit (dof_handler.n_dofs()) in assemble_system().

However, I was worried that these three values might not be substituted 
with new values(after every iteration).

So I added these three reinit in assemble_system().

The two results are a little bit different as expected.

When 3 reinit are not added...
   Number of active cells:   5120
   Number of degrees of freedom: 20609
  iteration 0
   ||u_{k}-u_{k-1}||42.5271
  iteration 1
   ||u_{k}-u_{k-1}||1.45458
  iteration 2
   ||u_{k}-u_{k-1}||0.899201
  iteration 3
   ||u_{k}-u_{k-1}||0.646361
  iteration 4
   ||u_{k}-u_{k-1}||0.48
  iteration 5
   ||u_{k}-u_{k-1}||0.404479
  iteration 6
   ||u_{k}-u_{k-1}||0.33739
  iteration 7
   ||u_{k}-u_{k-1}||0.287838
  iteration 8
   ||u_{k}-u_{k-1}||0.249859
  iteration 9
   ||u_{k}-u_{k-1}||0.219908

When added...

   Number of active cells:   5120
   Number of degrees of freedom: 20609
  iteration 0
   ||u_{k}-u_{k-1}||42.5271
  iteration 1
   ||u_{k}-u_{k-1}||5.2969
  iteration 2
   ||u_{k}-u_{k-1}||2.49065
  iteration 3
   ||u_{k}-u_{k-1}||1.45435
  iteration 4
   ||u_{k}-u_{k-1}||0.947504
  iteration 5
   ||u_{k}-u_{k-1}||0.66352
  iteration 6
   ||u_{k}-u_{k-1}||0.490408
  iteration 7
   ||u_{k}-u_{k-1}||0.378107
  iteration 8
   ||u_{k}-u_{k-1}||0.301494
  iteration 9
   ||u_{k}-u_{k-1}||0.246974

I'm confused that the error increases rather, as you can see.

I don't know what makes this difference and which one is right.

Thank you

Kyusik.

-- 
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] Isoparametric testing

[Lukas Bystricky]

I'm trying to test the convergence rates of isoparametric elements for the 
Poisson equation. To do this I have a circular domain and an exact 
solution. My understanding is that with linear mappings I should not get 
optimal convergence in the L2 and H1 norms for elements higher than degree 
1. However when I test this, I recover the optimal accuracy with quadratic, 
cubic and quartic elements. 

Here is my code to construct the mesh (using version 8.2.1)

Point<2> center (0,0);
double radius = 1;
GridGenerator::hyper_ball(
mesh, center, radius);
static const HyperBallBoundary boundary_description(center);
mesh.set_boundary(0, boundary_description);

I think that what's happening is that my computational domain is a polygon 
(using linear mappings) and when I apply the exact solution as a boundary 
condition deal.ii knows the exact value at all points along my polygonal 
boundary. In reality, it should only know the values on the circle and do 
some interpolation to figure out the values along the polygon. Here is my 
code to apply the boundary conditions:

Am I correct in this assumption? If so, is there a good way to test the 
convergence rates and demonstrate sub-optimal convergence for linear 
mappings?

-- 
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] Isoparametric testing

[Lukas Bystricky]

I'm trying to test the convergence rates of isoparametric elements for the 
Poisson equation. To do this I have a circular domain and an exact 
solution. My understanding is that with linear mappings I should not get 
optimal convergence in the L2 and H1 norms for elements higher than degree 
1. However when I test this, I recover the optimal accuracy with quadratic, 
cubic and quartic elements. 

Here is my code to construct the mesh (using version 8.2.1)

Point<2> center (0,0);
double radius = 1;
GridGenerator::hyper_ball(mesh, center, radius);
static const HyperBallBoundary boundary_description(center);
mesh.set_boundary(0, boundary_description);

I think that what's happening is that my computational domain is a polygon 
(using linear mappings) and when I apply the exact solution as a boundary 
condition deal.ii knows the exact value at all points along my polygonal 
boundary. In reality, it should only know the values on the circle and do 
some interpolation to figure out the values along the polygon. Here is my 
code to apply the boundary conditions:

std::map boundary_values;
VectorTools::interpolate_boundary_values (mapping, dof_handler, 0, 
*boundary_function, boundary_values);

MatrixTools::apply_boundary_values (boundary_values, stiffness_matrix,
  solution, rhs);

Am I correct in this assumption? If so, is there a good way to test the 
convergence rates and demonstrate sub-optimal convergence for lower order 
elements?



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

Re: [deal.II] bug in program only with adaptive mesh refinement

[Daniel Shapero]

>
> OK, so F_i would then be (f(u),\varphi_i), right? 
>
> I think the question is whether you do or do not apply 
> ConstraintMatrix::condense to F. I never really know whether you do or 
> do not have to do that when you want to use F that way, so I typically 
> just assemble the (scalar number) 
>(f(u),v) 
> directly. 
>

Hi Wolfgang, a brief update -- before, I was using 
ConstraintMatrix::distribute_local_to_global in the assembly of both the 
derivative vector and the Hessian matrix. I can see how that would make the 
computed vector no longer represent the derivative of a functional, since 
you're changing entries more to make the hanging node constraints come out 
correctly. To see if this was the issue, I instead changed the code to add 
up contributions to the global vector/matrix directly, and only reconcile 
the constraints in the nonlinear solver routine. All the other unit tests 
still pass, so this change didn't break the nonlinear solver for the 
unrefined case. Still, the errors in the local linear approximation to P 
aren't decreasing for the refined mesh.

I'll try your approach about assembling (f(u), v) directly and report back. 
If you have other suggestions, I'm all ears; this is a pretty serious 
roadblock for what I want to do next.

-- 
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: any suggestion how to run same code repeatedly with different boundary condition ?

[Jaekwang Kim]

Thank you !! solve the problem with the second method! 

I really appreciate this!

Jaekwang Kim 

2017년 3월 4일 토요일 오후 1시 37분 40초 UTC-6, Jean-Paul Pelteret 님의 말:
>
> Hi Jaekwang,
>
> The error message is pointing you to the fact that when you call run() a 
> second time you end up trying to rebuild the triangulation on top of the 
> one that you have the first time run() was called. There are a few ways 
> that you can work around this. Here are two suggestions:
>
> 1. Create a new instance of the laplace_problem each time you change the 
> boundary value
>
>   {
>
>bvalue = 1;
>
>Step4<2> laplace_problem_2d;
>
>laplace_problem_2d.run ();
>   }
>
>   {  
>
>bvalue = 2;
>   Step4<2> laplace_problem_2d;
>
>laplace_problem_2d.run ();
>
>  }
>
>
> 2. Separate the run() and make_grid() calls.
>
>   {
>
>Step4<2> laplace_problem_2d;
>   laplace_problem_2d.make_grid ();
>
>  
>bvalue = 1;
>laplace_problem_2d.run ();//Remove make_grid() from run()
>
>  
>
>bvalue = 2;
>
>laplace_problem_2d.run ();
>
>  
>
>  }
>
> You could of change your boundary value in other way, such as in a for 
> loop or using a MultipleParameterLoop 
> 
> .
>
> I hope this helps you.
>
> Regards,
> Jean-Paul
>
>

-- 
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: Question regarding P:D:Tria

[RAJAT ARORA]

Hello Daniel,

Thanks for the prompt reply. Yes, I thought about this as deal.ii uses 
P4est which does not support anisotropic refinement.

I wanted to confirm that this was indeed the case and I am not missing 
something. 


Thanks again for the clarification.

On Sunday, March 5, 2017 at 7:48:56 AM UTC-5, Daniel Arndt wrote:
>
> Rajat,
>
>
> I am using the GridGenerator::hyper_rectangle function to create a 
>> rectangular grid with 1000 X 500 X 1 (0.5 million) elements.
>>
>> My question is: - since the coarsest mesh contains 0.5 million elements, 
>> does that mean all the processors have 0.5 million elements?
>> (I think the answer is yes, I read that the deal.ii (p4est) stores the 
>> coarsest mesh on all processors.)
>>
> Yes, that is true. 
>
> I need a mesh (overall) which should be rectangular and should have 1000 X 
>> 500 X 1 elements. Note: It should only have 1 element in the z-direction.
>> Is there any way such that I can have this mesh by not replicating the 
>> whole mesh on all processors as it takes a lot of memory. I thought of 
>> making a 2D mesh and then refining it and then extruding it using 
>> extrude::triangulation, but extrude triangulation does not work on refined 
>> meshes.
>>
> As you already said, we (currently) require to store the same coarse mesh 
> on all processes and only isotropic refinement is allowed. This means that 
> you can get 1 element in z-direction if you don't refine. There is some 
> work [1] on a Triangulation that allows you to specify the distribution of 
> cells across processes yourself without the need of a common coarse mesh.
>
> Best,
> Daniel
>
> [1] https://github.com/dealii/dealii/pull/3956
>

-- 
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: Question regarding P:D:Tria

[Daniel Arndt]

Rajat,


I am using the GridGenerator::hyper_rectangle function to create a 
> rectangular grid with 1000 X 500 X 1 (0.5 million) elements.
>
> My question is: - since the coarsest mesh contains 0.5 million elements, 
> does that mean all the processors have 0.5 million elements?
> (I think the answer is yes, I read that the deal.ii (p4est) stores the 
> coarsest mesh on all processors.)
>
Yes, that is true. 

I need a mesh (overall) which should be rectangular and should have 1000 X 
> 500 X 1 elements. Note: It should only have 1 element in the z-direction.
> Is there any way such that I can have this mesh by not replicating the 
> whole mesh on all processors as it takes a lot of memory. I thought of 
> making a 2D mesh and then refining it and then extruding it using 
> extrude::triangulation, but extrude triangulation does not work on refined 
> meshes.
>
As you already said, we (currently) require to store the same coarse mesh 
on all processes and only isotropic refinement is allowed. This means that 
you can get 1 element in z-direction if you don't refine. There is some 
work [1] on a Triangulation that allows you to specify the distribution of 
cells across processes yourself without the need of a common coarse mesh.

Best,
Daniel

[1] https://github.com/dealii/dealii/pull/3956

-- 
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: Access specific element within a distributed triangulation

[Daniel Arndt]

Seyed,
 

> But is it necessary to #include  ?
> Since deal.II won't recognize functions such as MPI_Comm_rank 
> (MPI_COMM_WORLD, 
> &myrank) for instance or MPI_Reduce.
> Is it implemented in Utilities::MPI::something ? 
>
If you are using a parallel Triangulation you automatically include 
"". You can be sure that the comiler will tell you if you need to 
include that file explicitly.

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] can we bring variables from outside class...?

[Krzysztof Bzowski]

Hi,
the simplest solution would be just parametrize BoundaryValues 
constructor. It would look like:


template 
class BoundaryValues : public Function
{
double param;
public:

  BoundaryValues (double param) : Function() {
 this->param = param;
  }

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

};

Now you can use param inside any method in BoundaryValues class.

interpolate_boundary_values will look like:
VectorTools::interpolate_boundary_values (dof_handler, 0, 
BoundaryValues(*u_main), boundary_values);




I always learning a lot from you
For my problem, I have to impose different boundary values,
so for example if in step- 4, we use class derived from Function 
for

making boundary condition.

but.. is there any way that I can bring some values from mainclass? to 
be

more specific



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