Re: [deal.II] Poisson test code fails for a fe-degree >= 2

[Wolfgang Bangerth]

On 4/22/19 3:13 AM, 'Maxi Miller' via deal.II User Group wrote:
> How exactly can I test the examples after modifying them, other than by 
> running them?

In your case, the changes you want to make should not make a difference in the 
output of the program. (You're just replacing a fixed number by something like 
fe.degree+1, but these should be the same.) So if you run the program before 
and after the change and the output (on the screen, and to within eyeball norm 
also the graphical output) is the same, then that's a good enough test.

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.

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

['Maxi Miller' via deal.II User Group]

How exactly can I test the examples after modifying them, other than by 
running them?

Am Sonntag, 7. April 2019 23:56:49 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/7/19 1:54 PM, 'Maxi Miller' via deal.II User Group wrote: 
> > Could do that, but not until Easter (Deadlines before that), if that is 
> ok? 
>
> Any time is good for contributions! 
> Cheers 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

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

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

[Wolfgang Bangerth]

On 4/7/19 1:54 PM, 'Maxi Miller' via deal.II User Group wrote:
> Could do that, but not until Easter (Deadlines before that), if that is ok?

Any time is good for contributions!
Cheers
  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.

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

['Maxi Miller' via deal.II User Group]

Could do that, but not until Easter (Deadlines before that), if that is ok?

Am Freitag, 5. April 2019 17:39:30 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/5/19 9:34 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Oh, thanks! That happens if I start from an example, and don't write 
> > everything from scratch... 
>
> Since I stepped in that trap myself a couple of weeks ago and had to 
> have Timo rescue me, I'm thinking that we ought to just change this in 
> all tutorial programs where the quadrature formula degree is not tied to 
> the FE degree. Would you want to write up a pull request to make such a 
> change? 
>
> Best 
>   Wolfgang 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

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

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

[Wolfgang Bangerth]

On 4/5/19 9:34 AM, 'Maxi Miller' via deal.II User Group wrote:
> Oh, thanks! That happens if I start from an example, and don't write 
> everything from scratch...

Since I stepped in that trap myself a couple of weeks ago and had to 
have Timo rescue me, I'm thinking that we ought to just change this in 
all tutorial programs where the quadrature formula degree is not tied to 
the FE degree. Would you want to write up a pull request to make such a 
change?

Best
  Wolfgang

-- 

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.

Re: [deal.II] Poisson test code fails for a fe-degree >= 2

['Maxi Miller' via deal.II User Group]

Oh, thanks! That happens if I start from an example, and don't write 
everything from scratch...

Am Donnerstag, 4. April 2019 12:25:10 UTC+2 schrieb Martin Kronbichler:
>
> Dear Maxi,
>
> The problem is the quadrature formula you use in 'assemble_system' and 
> 'calculate_residual'. You should use 'QGauss 
> quadrature_formula(fe.degree+1)`. I am surprised that quadratic elements 
> even worked because the Poisson matrix for d>=2 should have some rank 
> deficiency with under-integration as far as I know.
>
> Best,
> Martin
> On 04.04.19 11:31, 'Maxi Miller' via deal.II User Group wrote:
>
> I wrote a test program, based on step-5, for testing unrelated things, and 
> checked the results. But when I increased the degree of the finite elements 
> from 1 to 3 (or higher), the results were wrong, while the change from 1 to 
> 2 just improved convergence, as expected. Did I make a mistake in my code 
> here?
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to 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] Poisson test code fails for a fe-degree >= 2

[Martin Kronbichler]

Dear Maxi,

The problem is the quadrature formula you use in 'assemble_system' and
'calculate_residual'. You should use 'QGauss
quadrature_formula(fe.degree+1)`. I am surprised that quadratic elements
even worked because the Poisson matrix for d>=2 should have some rank
deficiency with under-integration as far as I know.

Best,
Martin

On 04.04.19 11:31, 'Maxi Miller' via deal.II User Group wrote:
> I wrote a test program, based on step-5, for testing unrelated things,
> and checked the results. But when I increased the degree of the finite
> elements from 1 to 3 (or higher), the results were wrong, while the
> change from 1 to 2 just improved convergence, as expected. Did I make
> a mistake in my code here?
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to dealii+unsubscr...@googlegroups.com
> .
> For more options, visit https://groups.google.com/d/optout.

-- 
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] Poisson test code fails for a fe-degree >= 2

['Maxi Miller' via deal.II User Group]

I wrote a test program, based on step-5, for testing unrelated things, and 
checked the results. But when I increased the degree of the finite elements 
from 1 to 3 (or higher), the results were wrong, while the change from 1 to 
2 just improved convergence, as expected. Did I make a mistake in my code 
here?

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

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */


// @sect3{Include files}

// Again, the first few include files are already known, so we won't comment
// on them:
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

// This one is new. We want to read a triangulation from disk, and the class
// which does this is declared in the following file:
#include 

// We will use a circular domain, and the object describing the boundary of it
// comes from this file:
#include 

// This is C++ ...
#include 
#include 


// Finally, this has been discussed in previous tutorial programs before:
using namespace dealii;

template 
class Solution : public Function
{
public:
Solution() : Function(1)
{

}

virtual double value(const Point &p, const unsigned int component) const override;
virtual Tensor<1, dim> gradient(const Point &p, const unsigned int component) const override;

private:
};

template 
double Solution::value(const Point &p, const unsigned int) const
{
const double x = p[0];
const double y = p[1];
return sin(M_PI * x) * sin(M_PI * y);
}

template
Tensor<1, dim> Solution::gradient(const Point &p, const unsigned int) const
{
Tensor<1, dim> return_value;
AssertThrow(dim == 2, ExcNotImplemented());

const double x = p[0];
const double y = p[1];
return_value[0] = M_PI * cos(M_PI * x) * sin(M_PI * y);
return_value[1] = M_PI * cos(M_PI * y) * sin(M_PI * x);
return return_value;
}


// @sect3{The Step5 class template}

// The main class is mostly as in the previous example. The most visible
// change is that the function make_grid_and_dofs has been
// removed, since creating the grid is now done in the run
// function and the rest of its functionality is now in
// setup_system. Apart from this, everything is as before.
template 
class Step5
{
public:
Step5();
~Step5();
void run();

private:
void setup_system();
void assemble_system();
void calculate_residual_multi_val(const  Vector &epsilon_v, const Vector &u, Vector &residual);
void calculate_residual(const Vector &u, Vector &residual);
void solve();
void solve_with_matrix_jacobian();
void process_solution();
void output_results(const unsigned int cycle);

void assemble_rhs(double &return_value,
  const double factor,
  const double &val_u,
  const Tensor<1, dim> &grad_u,
  const Point p,
  const double &fe_values_value,
  const Tensor<1, dim> &fe_gradient_value,
  const double JxW_value);

Triangulation triangulation;
FE_Q  fe;
DoFHandlerdof_handler;

SparsityPattern  sparsity_pattern;
SparseMatrix system_matrix;

AffineConstraints hanging_node_constraints;

ConvergenceTable convergence_table;

Vector solution;
Vector system_rhs;
};


// @sect3{Working with nonconstant coefficients}

// In step-4, we showed how to use non-constant boundary values and right hand
// side.  In this example, we want to use a variable coefficient in the
// elliptic operator instead. Since we have a function which just depends on
// the point in space we can do things a bit more simply and use a plain
// function instead of inheriting from Function.

// This is the implementation of the coefficient function for a single
// p