[deal.II] Re: How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Dear Prof. Arndt,

I have used SolutionTransfer as you say, but if I set a phi_0 fixed, then 
project it into finite element space and get a vector phi_0_h, then I get 
phi_0_h/2, phi_0_h/4, phi_0_h/8 by using SolutionTransfer, but the norm 
of (phi_0_h- phi_0_h/2), (phi_0_h/2- phi_0_h/4), (phi_0_h/4- phi_0_h/8) are 
so big, and I cannot get exact value when I compute the convergence rate. 
So I wander how to deal with it? 

And the attached is the test code, and the results is:
phi_0_norm_58:  0.335377  phi_0_norm_68:  0.317158  phi_0_norm_78:  0.277344

Thank you very much!

Best,
Chucui

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include 

// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include 

// This is C++:
#include 
#include 
#include 
#include 
#include 
using namespace std;
// As in all programs, the namespace dealii is included:
namespace Step22
{
  using namespace dealii;


  template 
  class StokesProblem
  {
  public:
StokesProblem (const unsigned int phi_degree);
void run ();
  private:
void setup_dofs_pbc_rf (const int refine_number);
void phi_0_trans ();
 
double norm_compute_phi_0 (const Vector solu_1,
 const Vector solu_2); 
void create_mesh();
void refine_mesh ();

const unsigned int   degree;

Triangulation   triangulation_active;
Triangulation   triangulation;
FE_Qfe_phi;

DoFHandler  dof_handler_phi;
ConstraintMatrix constraints_phi;
Vector phi_0, phi_0_5, phi_0_6, phi_0_7, phi_0_8, tmp_phi_0_5, tmp_phi_0_6, tmp_phi_0_7, tmp_phi_0_8;
 
double time_step, time;
const double A, B, sita_0, eps, m, eps_4, tau, k_0, C_L, D_u, L_1, L_2, L_0, T_M, C_0, C_1, eps_F, e_L_T_M;//, s;
unsigned int timestep_number = 1;
  
  };

  
  template 
  StokesProblem::StokesProblem (const unsigned int phi_degree)
:
degree (phi_degree),
fe_phi (degree),
dof_handler_phi (triangulation),
time_step (1e-4),
timestep_number (1),
A (100.0), 
B (1.0), 
sita_0 (1*numbers::PI/4), 
eps (0.01), 
m (4.0), 
eps_4 (1./15.0), 
tau (3.0), 
k_0 (0.0), 
C_L (0.6), 
D_u (0.0001), 
L_1 (0.0), 
L_2 (0.0), 
L_0 (0.5), 
T_M (1.0), 
C_0 (1.0), 
C_1 (1.0),
eps_F (1.0),
e_L_T_M (0.0)
  {}  


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

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

  };
  
  
  template 
  double InitialValues_phi::value (const Point   &p,
   const unsigned int component) const
  {
  const double C_0 = 1.0,
   C_1 = 1.0,
   A = 100.0;
  const double delta = 0.1, r0 = std::sqrt(0.018);
  double phi_value = 0, e_value = 0, q_value = 0, s_value = 0, r = 0, T_value = 0;
   r = std::sqrt((p[0] - 0.28125) * (p[0] - 0.28125) + (p[1] - 0.28125) * (p[1] - 0.28125)) ;
 phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
 //phi_value = 0.5 * std::cos(p[0]) * std::cos(p[1]) +0.5;
 return phi_value;
  
  }  

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

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

  };


  template 
  double
  InitialValuesConv::value (const Point  &p,
 const unsigned int component) const

[deal.II] Re: How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Dear Prof. Arndt,

I have used SolutionTransfer as you say, but if I set a phi_0 fixed, then 
project it into finite element space and get a vector phi_0_h, then I get 
phi_0_h/2, phi_0_h/4, phi_0_h/8 by using SolutionTransfer, but the norm 
of (phi_0_h- phi_0_h/2), (phi_0_h/2- phi_0_h/4), (phi_0_h/4- phi_0_h/8) are 
so big, and I cannot get exact value when I compute the convergence rate. 
So I wander how to deal with it? 

And the attached is the test code, and the results is:
phi_0_norm_58:  0.335377  phi_0_norm_68:  0.317158  phi_0_norm_78:  0.277344

Thank you very much!

Best,
Chucui

-- 
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: How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Dear Prof.Arndt,

Thank you very much! I understand what you say. It helps me a lot!

Best,
Chucui

在 2019年2月21日星期四 UTC+8下午7:31:20，Daniel Arndt写道：
>
> Chucui,
>
>
> For Question 1, I write a code to compute the L2 norm of (solution_1 - 
>> solution_2):
>> [...]
>>
>> Is that right?
>>
> You can just take the difference of the two variables and use 
> VectorTools::integrate_difference with a Functions::ZeroFunction as "exact 
> solution".
>  
>  
>
>>
>> For question 2, as the 2 solution vectors have different sizes, we don't 
>> have the same cells, so the cell loop cannot compute altogether:
>> typename DoFHandler::active_cell_iterator
>> cell_1 = dof_handler_1.begin_active(),
>> endc_1 = dof_handler_1.end(); 
>> typename DoFHandler::active_cell_iterator
>> cell_2 = dof_handler_2.begin_active();  
>> for (; cell_1!=endc_2; ++cell_1, ++cell_2)
>>   {}
>>
>> I wander how to solve this problem?
>>
> You need to interpolate the solution computed on the coarser mesh to the 
> finer mesh using SolutionTransfer(
> https://www.dealii.org/current/doxygen/deal.II/classSolutionTransfer.html, 
> explained e.g. in 
> https://www.dealii.org/current/doxygen/deal.II/step_26.html). Then, you 
> can take the difference just as in your first question.
>
> 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] How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Dear Jean-Paul,

Thank you very much for your quick reply! I will read it now! 

Best,
Chucui

在 2019年2月21日星期四 UTC+8下午5:17:17，Jean-Paul Pelteret写道：
>
> Hi Chucui,
>
> I can only offer you a very brief reply right now: You might want to look 
> into the ConvergenceTable 
>  
> class 
> to help you with (1). The introduction to the class directs you to some 
> tutorials that indicate its use.
>
> Best,
> Jean-Paul
>
> On 21 Feb 2019, at 10:11, chucu...@gmail.com  wrote:
>
> Hi all,
>
> I have 2 questions of computing convergence rate of L2 norm of error:
>
> 1. If I don't have exact solution, does deal.ii have functions to compute 
> convergence rate directly? (Like VectorTools::integrate_difference, which 
> deals with the problem with exact solution  )  
>   
> 2. If  I want to compute convergence rate of error of solutions with 
> different globally refinement (u_{h}-u_{h/2}, u_{h/2}-u_{h/4}, and so on ), 
> does deal.ii have  functions to compute convergence rate directly?
> if not, how to deal with computing u_{h}-u_{h/2}, where the 2 numerical 
> solution have different size of vector?
>
> Thanks in advance!
>
> Best,
> Chucui
>
> -- 
> 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.

[deal.II] Re: How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Hi, all,

If there is no function to compute convergence above, I need to write the 
code by myself.

For Question 1, I write a code to compute the L2 norm of (solution_1 - 
solution_2):

template 
  double StokesProblem::norm_compute (const Vector solu_1,
  const Vector solu_2)
  {
QGauss   quadrature_formula(degree+2);

FEValues fe_values (fe, quadrature_formula,
 update_values|
 update_quadrature_points  |
 update_JxW_values |
 update_gradients);
  
const unsigned int   dofs_per_cell   = fe.dofs_per_cell;

const unsigned int   n_q_points  = quadrature_formula.size();
FullMatrix   local_con_rate_matrix (dofs_per_cell, 
dofs_per_cell);
std::vector local_dof_indices (dofs_per_cell); 

std::vector solu_1_values(n_q_points);
std::vector solu_2_values(n_q_points);

double norm_12 = 0, norm_24 = 1,norm_12_square = 0, norm_24_square = 1, 
rate_124 = 0, con_rate = 0; 
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();   
for (; cell!=endc; ++cell)
  {
fe_values.reinit (celli);  
fe_values.get_function_values (solu_1, solu_1_values);
fe_values.get_function_values (solu_2, solu_2_values);
for (unsigned int q=0; q::active_cell_iterator
cell_1 = dof_handler_1.begin_active(),
endc_1 = dof_handler_1.end(); 
typename DoFHandler::active_cell_iterator
cell_2 = dof_handler_2.begin_active();  
for (; cell_1!=endc_2; ++cell_1, ++cell_2)
  {}

I wander how to solve this problem?

Thanks in advance!
Best,
Chucui

-- 
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] How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

[chucui1016]

Hi all,

I have 2 questions of computing convergence rate of L2 norm of error:

1. If I don't have exact solution, does deal.ii have functions to compute 
convergence rate directly? (Like VectorTools::integrate_difference, which 
deals with the problem with exact solution  )  
  
2. If  I want to compute convergence rate of error of solutions with 
different globally refinement (u_{h}-u_{h/2}, u_{h/2}-u_{h/4}, and so on ), 
does deal.ii have  functions to compute convergence rate directly?
if not, how to deal with computing u_{h}-u_{h/2}, where the 2 numerical 
solution have different size of vector?

Thanks in advance!

Best,
Chucui

-- 
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] different results when compute integration with adaptive mesh

[chucui1016]

Dear Jean-Paul,

Thank you very much for you answer! It helps me a lot. The attached is the 
implemention of your idea, and the result is:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
u_square:  0.0305621  grad_u_square:  0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
u_square:  0.0425496  grad_u_square:  0.315907
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
u_square:  0.0489012  grad_u_square:  0.341164
Cycle 3:
   Number of active cells:   200
   Number of degrees of freedom: 921
u_square:  0.0521656  grad_u_square:  0.353416
Cycle 4:
   Number of active cells:   440
   Number of degrees of freedom: 2017
u_square:  0.0540783  grad_u_square:  0.36169
Cycle 5:
   Number of active cells:   956
   Number of degrees of freedom: 4425
u_square:  0.0544531  grad_u_square:  0.363144
Cycle 6:
   Number of active cells:   1916
   Number of degrees of freedom: 8993
u_square:  0.0552556  grad_u_square:  0.366492
Cycle 7:
   Number of active cells:   3860
   Number of degrees of freedom: 18353
u_square:  0.0555461  grad_u_square:  0.36768


Thank you very much! And also thank Prof.Bangerth and Prof.Arndt! And I 
will think about the constraints much more.

Best,
Chucui 


在 2019年2月12日星期二 UTC+8下午4:01:10，chucu...@gmail.com写道：
>
> Dear Jean-Paul,
>
> Thank you very much for your detailed answer! I will write it as you say 
> right now and report the results as soon as possible !
>
> Best,
> Chucui
>
> 在 2019年2月12日星期二 UTC+8下午3:26:34，Jean-Paul Pelteret写道：
>>
>> Dear Chucui, 
>>
>> Too add to all of the other comments, you can "sanity check” the result 
>> that you’re getting using this approach by performing the same calculations 
>> another way and comparing results. Specifically, you could easily calculate 
>> these norms by looping over all cells and integrating the quantities 
>> yourself. As pseudocode it would look like this: 
>> 1. Loop over all cells in DoFHandler 
>> 2. Reinit FEValues with the update_values and update_gradients flags 
>> enabled 
>> 3. Compute local solution and its gradient at the cell quadrature points 
>> using fe_values.get_solution_[values/gradients]() 
>> 4. Loop over all quadrature points 
>> 5. u_square +=  local_solution[q_point]*local_solution[q_point] * 
>> fe_values.JxW() ; grad_u_square +=  local_gradient[q_point]* 
>> local_gradient[q_point] * fe_values.JxW() 
>>
>> Maybe this approach could help you debug your issue. 
>> Best, 
>> Jean-Paul 
>>
>> > On 12 Feb 2019, at 06:49, Wolfgang Bangerth  
>> wrote: 
>> > 
>> > On 2/11/19 12:20 AM, chucu...@gmail.com wrote: 
>> >> Dear Prof. Bangerth: 
>> >> 
>> >> Thank you very much for your quick reply! 
>> >> 
>> >> When I apply hanging node constraints to my matrix, as the step-6 
>> >> says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
>> >> I make 4 steps: 
>> >>   1.Create a Constraints Class: 
>> >> | 
>> >> ConstraintMatrixconstraints; 
>> >> | 
>> >> 
>> >>   2.Fill this object using the function 
>> >> DoFTools::make_hanging_node_constraints() to ensure continuity of the 
>> elements 
>> >> of the finite element space. 
>> >> | 
>> >> DoFTools::make_hanging_node_constraints (dof_handler, 
>> >>  constraints); 
>> >> | 
>> >> 
>> >>   3.Copy the local contributions to the matrix and right hand side 
>> into the 
>> >> global objects: 
>> >> | 
>> >> constraints.distribute_local_to_global (cell_matrix, 
>> >> cell_rhs, 
>> >> local_dof_indices, 
>> >> system_matrix, 
>> >> system_rhs); 
>> >> constraints.distribute_local_to_global (cell_volume_matrix, 
>> >> local_dof_indices, 
>> >> volume_matrix); 
>> >> constraints.distribute_local_to_global (cell_gradient_matrix, 
>> >> local_dof_indices, 
>> >> gradient_matrix); 
>> >> | 
>> >> 
>> >> 
>> >>   4.Make sure that the degrees of "freedom" located on hanging nodes 
>> get 
>> >> their correct (constrained) value: 
>> >> | 
>> >> constraints.distribute (solution); 
>> >> | 
>> > 
>> > Constraints are funny and sometimes require deep thought about what 
>> exactly 
>> > they mean. What happens if you don't apply constraints to the 
>> > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
>> elements 1:1 
>> > into the matrix, without the 'constraints' object? (Like in step-4.) 
>> > 
>> > Alternatively, you could do as you show above and after calling 
>> > 'constraint.distribute(solution)' you manually set all constrained 
>> degrees of 
>> > freedom in the solution vector to zero. That's not what you want to use 
>> f

Re: [deal.II] different results when compute integration with adaptive mesh

[chucui1016]

Dear Jean-Paul,

Thank you very much for your detailed answer! I will write it as you say 
right now and report the results as soon as possible !

Best,
Chucui

在 2019年2月12日星期二 UTC+8下午3:26:34，Jean-Paul Pelteret写道：
>
> Dear Chucui, 
>
> Too add to all of the other comments, you can "sanity check” the result 
> that you’re getting using this approach by performing the same calculations 
> another way and comparing results. Specifically, you could easily calculate 
> these norms by looping over all cells and integrating the quantities 
> yourself. As pseudocode it would look like this: 
> 1. Loop over all cells in DoFHandler 
> 2. Reinit FEValues with the update_values and update_gradients flags 
> enabled 
> 3. Compute local solution and its gradient at the cell quadrature points 
> using fe_values.get_solution_[values/gradients]() 
> 4. Loop over all quadrature points 
> 5. u_square +=  local_solution[q_point]*local_solution[q_point] * 
> fe_values.JxW() ; grad_u_square +=  local_gradient[q_point]* 
> local_gradient[q_point] * fe_values.JxW() 
>
> Maybe this approach could help you debug your issue. 
> Best, 
> Jean-Paul 
>
> > On 12 Feb 2019, at 06:49, Wolfgang Bangerth  > wrote: 
> > 
> > On 2/11/19 12:20 AM, chucu...@gmail.com  wrote: 
> >> Dear Prof. Bangerth: 
> >> 
> >> Thank you very much for your quick reply! 
> >> 
> >> When I apply hanging node constraints to my matrix, as the step-6 
> >> says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
> >> I make 4 steps: 
> >>   1.Create a Constraints Class: 
> >> | 
> >> ConstraintMatrixconstraints; 
> >> | 
> >> 
> >>   2.Fill this object using the function 
> >> DoFTools::make_hanging_node_constraints() to ensure continuity of the 
> elements 
> >> of the finite element space. 
> >> | 
> >> DoFTools::make_hanging_node_constraints (dof_handler, 
> >>  constraints); 
> >> | 
> >> 
> >>   3.Copy the local contributions to the matrix and right hand side into 
> the 
> >> global objects: 
> >> | 
> >> constraints.distribute_local_to_global (cell_matrix, 
> >> cell_rhs, 
> >> local_dof_indices, 
> >> system_matrix, 
> >> system_rhs); 
> >> constraints.distribute_local_to_global (cell_volume_matrix, 
> >> local_dof_indices, 
> >> volume_matrix); 
> >> constraints.distribute_local_to_global (cell_gradient_matrix, 
> >> local_dof_indices, 
> >> gradient_matrix); 
> >> | 
> >> 
> >> 
> >>   4.Make sure that the degrees of "freedom" located on hanging nodes 
> get 
> >> their correct (constrained) value: 
> >> | 
> >> constraints.distribute (solution); 
> >> | 
> > 
> > Constraints are funny and sometimes require deep thought about what 
> exactly 
> > they mean. What happens if you don't apply constraints to the 
> > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
> elements 1:1 
> > into the matrix, without the 'constraints' object? (Like in step-4.) 
> > 
> > Alternatively, you could do as you show above and after calling 
> > 'constraint.distribute(solution)' you manually set all constrained 
> degrees of 
> > freedom in the solution vector to zero. That's not what you want to use 
> for 
> > visualization, error estimation, ... but it might work for the operation 
> > you're trying here. 
> > 
> > I'm pretty sure that both of these solutions should work individually 
> (but not 
> > in combination). But explaining why this is so would take a page or two, 
> I'm 
> > afraid. 
> > 
> > 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+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

Re: [deal.II] Re: different results when compute integration with adaptive mesh

[chucui1016]

The attached are codes:
int-test-1st-result.cc—— for the first result
int-test-2nd-result.cc—— for the second result

在 2019年2月12日星期二 UTC+8下午3:40:57，chucu...@gmail.com写道：
>
> Dear Prof. Bangerth,
>
> Thanks very much for your quick answer!
>
> > Constraints are funny and sometimes require deep thought about what 
> exactly 
> > they mean. What happens if you don't apply constraints to the 
> > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
> elements 1:1 
> > into the matrix, without the 'constraints' object? (Like in step-4.) 
> I rewrite my code as you say, and the results are:
> Cycle 0:
>Number of active cells:   20
>Number of degrees of freedom: 89
>  u square 0: 0.0305621  grad u square 0: 0.256542
> Cycle 1:
>Number of active cells:   44
>Number of degrees of freedom: 209
>  u square 1: 0.144303  grad u square 1: 0.49794
> Cycle 2:
>Number of active cells:   92
>Number of degrees of freedom: 449
>  u square 2: 0.101366  grad u square 2: 0.477137
> Cycle 3:
>Number of active cells:   188
>Number of degrees of freedom: 881
>  u square 3: 0.13305  grad u square 3: 0.529805
> Cycle 4:
>Number of active cells:   368
>Number of degrees of freedom: 1737
>  u square 4: 0.145489  grad u square 4: 0.565953
> Cycle 5:
>Number of active cells:   737
>Number of degrees of freedom: 3409
>  u square 5: 0.141771  grad u square 5: 0.574863
> Cycle 6:
>Number of active cells:   1436
>Number of degrees of freedom: 6705
>  u square 6: 0.171807  grad u square 6: 0.642486
> Cycle 7:
>Number of active cells:   2729
>Number of degrees of freedom: 12521
>  u square 7: 0.160819  grad u square 7: 0.625337
>
> The attached is the code ,which doesn't use 
> "constraints.distribute_local_to_global(...)" and 
> "DoFTools::make_hanging_node_constraints (dof_handler,  constraints);"
> if "DoFTools::make_hanging_node_constraints (dof_handler,  constraints);" 
> remain, the results:
> Cycle 0:
>Number of active cells:   20
>Number of degrees of freedom: 89
>  u square 0: 0.0305621  grad u square 0: 0.256542
> Cycle 1:
>Number of active cells:   44
>Number of degrees of freedom: 209
>  u square 1: 0.0822159  grad u square 1: 2.97512
> Cycle 2:
>Number of active cells:   92
>Number of degrees of freedom: 449
>  u square 2: 0.000404454  grad u square 2: 0.0534184
> Cycle 3:
>Number of active cells:   188
>Number of degrees of freedom: 937
>  u square 3: 9.40228e-05  grad u square 3: 0.0378517
> Cycle 4:
>Number of active cells:   368
>Number of degrees of freedom: 1793
>  u square 4: 7.66868e-05  grad u square 4: 0.0403114
> Cycle 5:
>Number of active cells:   704
>Number of degrees of freedom: 3321
>  u square 5: 5.9193e-05  grad u square 5: 0.0366146
> Cycle 6:
>Number of active cells:   1412
>Number of degrees of freedom: 6657
>  u square 6: 6.47763e-05  grad u square 6: 0.0455655
> Cycle 7:
>Number of active cells:   2696
>Number of degrees of freedom: 12477
>  u square 7: 2.50787e-05  grad u square 7: 0.0320763
>
>
>
> So I wander if it approximate the "grad u_exact square" and  "u_exact 
> square" correctly?  
> Because in my numerical experiments, the  "grad u_h square" and  "u_h 
> square" of numerical solution u_h are what we take care of. 
>
> So if the computions of "grad u_h square" and  "u_h square" are divided 
> from other work like visualization, error estimation, are these results 
> correct if change the using of constraints?
> And without "constraints..distribute_local_to_global(...)", the 
> information of boundary may not use correctly (for example, My boundary 
> condition is periodic boundary condition).
>
> Thank you very much!
>
> Best,
> Chucui
>
> 在 2019年2月12日星期二 UTC+8下午1:49:15，Wolfgang Bangerth写道：
>>
>> On 2/11/19 12:20 AM, chucu...@gmail.com wrote: 
>> > Dear Prof. Bangerth: 
>> > 
>> > Thank you very much for your quick reply! 
>> > 
>> > When I apply hanging node constraints to my matrix, as the step-6 
>> > says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
>> > I make 4 steps: 
>> >  1.Create a Constraints Class: 
>> > | 
>> > ConstraintMatrixconstraints; 
>> > | 
>> > 
>> >  2.Fill this object using the function 
>> > DoFTools::make_hanging_node_constraints() to ensure continuity of the 
>> elements 
>> > of the finite element space. 
>> > | 
>> > DoFTools::make_hanging_node_constraints (dof_handler, 
>> > constraints); 
>> > | 
>> > 
>> >  3.Copy the local contributions to the matrix and right hand side 
>> into the 
>> > global objects: 
>> > | 
>> >constraints.distribute_local_to_global (cell_matrix, 
>> >cell_rhs, 
>> >local_dof_indices, 
>> >system_matrix, 
>> >

Re: [deal.II] Re: different results when compute integration with adaptive mesh

[chucui1016]

Dear Prof. Bangerth,

Thanks very much for your quick answer!

> Constraints are funny and sometimes require deep thought about what 
exactly 
> they mean. What happens if you don't apply constraints to the 
> cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
elements 1:1 
> into the matrix, without the 'constraints' object? (Like in step-4.) 
I rewrite my code as you say, and the results are:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
 u square 0: 0.0305621  grad u square 0: 0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
 u square 1: 0.144303  grad u square 1: 0.49794
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
 u square 2: 0.101366  grad u square 2: 0.477137
Cycle 3:
   Number of active cells:   188
   Number of degrees of freedom: 881
 u square 3: 0.13305  grad u square 3: 0.529805
Cycle 4:
   Number of active cells:   368
   Number of degrees of freedom: 1737
 u square 4: 0.145489  grad u square 4: 0.565953
Cycle 5:
   Number of active cells:   737
   Number of degrees of freedom: 3409
 u square 5: 0.141771  grad u square 5: 0.574863
Cycle 6:
   Number of active cells:   1436
   Number of degrees of freedom: 6705
 u square 6: 0.171807  grad u square 6: 0.642486
Cycle 7:
   Number of active cells:   2729
   Number of degrees of freedom: 12521
 u square 7: 0.160819  grad u square 7: 0.625337

The attached is the code ,which doesn't use 
"constraints.distribute_local_to_global(...)" and 
"DoFTools::make_hanging_node_constraints (dof_handler,  constraints);"
if "DoFTools::make_hanging_node_constraints (dof_handler,  constraints);" 
remain, the results:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
 u square 0: 0.0305621  grad u square 0: 0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
 u square 1: 0.0822159  grad u square 1: 2.97512
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
 u square 2: 0.000404454  grad u square 2: 0.0534184
Cycle 3:
   Number of active cells:   188
   Number of degrees of freedom: 937
 u square 3: 9.40228e-05  grad u square 3: 0.0378517
Cycle 4:
   Number of active cells:   368
   Number of degrees of freedom: 1793
 u square 4: 7.66868e-05  grad u square 4: 0.0403114
Cycle 5:
   Number of active cells:   704
   Number of degrees of freedom: 3321
 u square 5: 5.9193e-05  grad u square 5: 0.0366146
Cycle 6:
   Number of active cells:   1412
   Number of degrees of freedom: 6657
 u square 6: 6.47763e-05  grad u square 6: 0.0455655
Cycle 7:
   Number of active cells:   2696
   Number of degrees of freedom: 12477
 u square 7: 2.50787e-05  grad u square 7: 0.0320763



So I wander if it approximate the "grad u_exact square" and  "u_exact 
square" correctly?  
Because in my numerical experiments, the  "grad u_h square" and  "u_h 
square" of numerical solution u_h are what we take care of. 

So if the computions of "grad u_h square" and  "u_h square" are divided 
from other work like visualization, error estimation, are these results 
correct if change the using of constraints?
And without "constraints..distribute_local_to_global(...)", the information 
of boundary may not use correctly (for example, My boundary condition is 
periodic boundary condition).

Thank you very much!

Best,
Chucui

在 2019年2月12日星期二 UTC+8下午1:49:15，Wolfgang Bangerth写道：
>
> On 2/11/19 12:20 AM, chucu...@gmail.com  wrote: 
> > Dear Prof. Bangerth: 
> > 
> > Thank you very much for your quick reply! 
> > 
> > When I apply hanging node constraints to my matrix, as the step-6 
> > says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
> > I make 4 steps: 
> >  1.Create a Constraints Class: 
> > | 
> > ConstraintMatrixconstraints; 
> > | 
> > 
> >  2.Fill this object using the function 
> > DoFTools::make_hanging_node_constraints() to ensure continuity of the 
> elements 
> > of the finite element space. 
> > | 
> > DoFTools::make_hanging_node_constraints (dof_handler, 
> > constraints); 
> > | 
> > 
> >  3.Copy the local contributions to the matrix and right hand side 
> into the 
> > global objects: 
> > | 
> >constraints.distribute_local_to_global (cell_matrix, 
> >cell_rhs, 
> >local_dof_indices, 
> >system_matrix, 
> >system_rhs); 
> >constraints.distribute_local_to_global (cell_volume_matrix, 
> >local_dof_indices, 
> >volume_matrix); 
> >constraints.distribute_local_to_global (cell_gradient_matrix, 
> >loca

Re: [deal.II] Re: different results when compute integration with adaptive mesh

[chucui1016]

Dear Prof Arndt,

Thanks for your quick answer!

> in the end you program is very similar to step-6. Do you observe the same 
behavior there as well?
Yes, my code add several lines below to step-6.cc directly:
...
  SparseMatrix volume_matrix;
  SparseMatrix gradient_matrix;
...
  volume_matrix.reinit (sparsity_pattern);
  gradient_matrix.reinit (sparsity_pattern);
...

 cell_volume_matrix(i,j) += (
 fe_values.shape_value(i,q_index) *
 fe_values.shape_value(j,q_index) *
 fe_values.JxW(q_index));
 cell_gradient_matrix(i,j) += (
 fe_values.shape_grad(i,q_index) *
 fe_values.shape_grad(j,q_index) *
 fe_values.JxW(q_index));
...
  constraints.distribute_local_to_global (cell_volume_matrix,
  local_dof_indices,
  volume_matrix);
  constraints.distribute_local_to_global (cell_gradient_matrix,
  local_dof_indices,
  gradient_matrix); 
...
  double int_test_volume = volume_matrix.matrix_norm_square(solution);
  double int_test_gradient = 
gradient_matrix.matrix_norm_square(solution);
  std::cout << " u square " << cycle <<": " << int_test_volume << "  
grad u square " << cycle <<": " << int_test_gradient  << std::endl;
and the results:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
 u square 0: 0.0305621  grad u square 0: 0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
 u square 1: 0.0438936  grad u square 1: 1.31255
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
 u square 2: 0.0505902  grad u square 2: 3.584
Cycle 3:
   Number of active cells:   200
   Number of degrees of freedom: 921
 u square 3: 0.0536183  grad u square 3: 6.60495
Cycle 4:
   Number of active cells:   440
   Number of degrees of freedom: 2017
 u square 4: 0.0555872  grad u square 4: 19.7488
Cycle 5:
   Number of active cells:   956
   Number of degrees of freedom: 4425
 u square 5: 0.0562045  grad u square 5: 48.0924
Cycle 6:
   Number of active cells:   1916
   Number of degrees of freedom: 8993
 u square 6: 0.0572988  grad u square 6: 117.562
Cycle 7:
   Number of active cells:   3860
   Number of degrees of freedom: 18353
 u square 7: 0.0577537  grad u square 7: 258.867


Thank you very much! And attached is the code of step-6.cc added lines 
above.

Best,
Chucui


在 2019年2月11日星期一 UTC+8下午5:52:56，Daniel Arndt写道：
>
> Chucui,
>
> in the end you program is very similar to step-6. Do you observe the same 
> behavior there as well?
> If not, try to find the changes that are responsible for the quickly 
> increasing norm.
> Do the results of your computation look reasonable?
>
> And the version of deal.ii I use is 8.5.1, which doesn't have the updated 
>> class AffineConstraints , but have Constraint_matrix Class. I don't 
>> know if this is the key point. Need I update my deal.ii?
>>
> No, that should not be related at all.
>
> 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.
/* -
 *
 * Copyright (C) 2000 - 2016 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 at
 * the top level of the deal.II distribution.
 *
 * -

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


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

// From the following include file we will import the declaration of
// H1-conforming finite element shape functions. This family of finite
/

Re: [deal.II] Re: different results when compute integration with adaptive mesh

[chucui1016]

Dear Prof. Bangerth:

Thank you very much for your quick reply! 

When I apply hanging node constraints to my matrix, as the step-6 
says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html
I make 4 steps:
1.Create a Constraints Class:
ConstraintMatrix constraints;

2.Fill this object using the function 
DoFTools::make_hanging_node_constraints() to ensure continuity of the 
elements of the finite element space.
DoFTools::make_hanging_node_constraints (dof_handler,
   constraints);

3.Copy the local contributions to the matrix and right hand side into 
the global objects:
  constraints.distribute_local_to_global (cell_matrix,
  cell_rhs,
  local_dof_indices,
  system_matrix,
  system_rhs);
  constraints.distribute_local_to_global (cell_volume_matrix,
  local_dof_indices,
  volume_matrix);
  constraints.distribute_local_to_global (cell_gradient_matrix,
  local_dof_indices,
  gradient_matrix); 


4.Make sure that the degrees of "freedom" located on hanging nodes get 
their correct (constrained) value:
 constraints.distribute (solution);

So I don't know which step I use incorrectly? And the version of deal.ii I 
use is 8.5.1, which doesn't have the updated class AffineConstraints 
, but have Constraint_matrix Class. I don't know if this is the key 
point. Need I update my deal.ii?

Thank you very much!

Best,
Chucui

在 2019年2月11日星期一 UTC+8下午1:40:07，Wolfgang Bangerth写道：
>
> On 2/10/19 8:20 PM, chucu...@gmail.com  wrote: 
> > 
> > And the results without adaptive mesh seem like approximate the exact 
> results 
> > correctly. So I am confused why the results of "grad u square" with 
> adaptive 
> > mesh are so strange? 
>
> I suspect that you forget to apply hanging node constraints to your 
> matrix, or 
> apply them incorrectly -- which is exactly why I asked you to try globally 
> refined meshes: they don't have hanging node constraints :-) 
>
> 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] Re: different results when compute integration with adaptive mesh

[chucui1016]

The attach is the code for refine global.

在 2019年2月11日星期一 UTC+8上午11:21:00，chucu...@gmail.com写道：
>
> Dear Prof. Bangerth,
>>
>
> Thanks for your quick reply! 
>
> >  You mean you are wondering why the "grad u square" term grows with the 
> size of the problem but "u square" does not? 
>
> Yes, that's what I am confused. Our u_exact have exact "grad u_exact 
> square" and  "u_exact square". And our numerical solution "u" and its  
> "grad u square" and  "u square" should be approximation of the exact one 
> though with adaptive mesh. So I think  "grad u square" should not have so 
> huge jump with different adaptive meshes, because they are all the 
> approximation of   "grad u_exact square".
>
> > I don't know either (and the one place I looked at in your code seems 
> correct to me), but what happens if you do global refinement? Does it 
> increase with a fixed ratio from one step to the next? 
>
> I change the code in attaching, and the output is:
> Cycle 0:
>Number of active cells:   20
>Number of degrees of freedom: 89
>  u square 0: 0.0305621  grad u square 0: 0.256542
> Cycle 1:
>Number of active cells:   80
>Number of degrees of freedom: 337
>  u square 1: 0.0489011  grad u square 1: 0.341164
> Cycle 2:
>Number of active cells:   320
>Number of degrees of freedom: 1313
>  u square 2: 0.0540201  grad u square 2: 0.361384
> Cycle 3:
>Number of active cells:   1280
>Number of degrees of freedom: 5185
>  u square 3: 0.0549874  grad u square 3: 0.365557
> Cycle 4:
>Number of active cells:   5120
>Number of degrees of freedom: 20609
>  u square 4: 0.0553762  grad u square 4: 0.367069
> Cycle 5:
>Number of active cells:   20480
>Number of degrees of freedom: 82177
>  u square 5: 0.0555421  grad u square 5: 0.367712
> Cycle 6:
>Number of active cells:   81920
>Number of degrees of freedom: 328193
>  u square 6: 0.0556186  grad u square 6: 0.367986
> Cycle 7:
>Number of active cells:   327680
>Number of degrees of freedom: 1311745
>  u square 7: 0.0556576  grad u square 7: 0.368119
>
>
> and the ratio 
>
> [image: ask-8-1.JPG]
> are:
> u_square:1.84102.40381.31491.22871.1168
> 0.9720
> grad_u_square:   2.06522.27661.46461.23361.28430.9891
>
> And the results without adaptive mesh seem like approximate the exact 
> results correctly. So I am confused why the results of "grad u square" with 
> adaptive mesh are so strange? Maybe the reason is that we refinre mesh 
> where the gradient of u is huge, so only "grad_u_square" can be very 
> different with changing adaptive meshes? If so, how to deal with this 
> problem? How to approximate   "grad u_ecaxt square" correctly with changing 
> adaptive mesh? Maybe change the value smaller of scale of refine cells and 
> coarsen cells, to make the change smaller?
>
> Thank you very much!
>
> Best,
> Chucui
>  
>
>

-- 
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) 2000 - 2016 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 at
 * the top level of the deal.II distribution.
 *
 * -

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


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

#include 

#include 

#include 

// Finally, this is as in previous programs:
using namespace dealii;

template 
class Step6
{
public:
  Step6 ();
  ~Step6 ();

  void run ();

private:
  void setup_system ();
  void assemble_system ();
  void solve ();
  void refine_grid ();
  void output_results (const unsigned int cycle) const;

  Triangulation   triangulation;

  DoFHandler  dof_handler;
  FE_Qfe;

  ConstraintMatrix constraints;

  SparsityPattern  sparsity_pattern;
  

[deal.II] Re: different results when compute integration with adaptive mesh

[chucui1016]

>
> Dear Prof. Bangerth,
>

Thanks for your quick reply! 

>  You mean you are wondering why the "grad u square" term grows with the 
size of the problem but "u square" does not? 

Yes, that's what I am confused. Our u_exact have exact "grad u_exact 
square" and  "u_exact square". And our numerical solution "u" and its  
"grad u square" and  "u square" should be approximation of the exact one 
though with adaptive mesh. So I think  "grad u square" should not have so 
huge jump with different adaptive meshes, because they are all the 
approximation of   "grad u_exact square".

> I don't know either (and the one place I looked at in your code seems 
correct to me), but what happens if you do global refinement? Does it 
increase with a fixed ratio from one step to the next? 

I change the code in attaching, and the output is:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
 u square 0: 0.0305621  grad u square 0: 0.256542
Cycle 1:
   Number of active cells:   80
   Number of degrees of freedom: 337
 u square 1: 0.0489011  grad u square 1: 0.341164
Cycle 2:
   Number of active cells:   320
   Number of degrees of freedom: 1313
 u square 2: 0.0540201  grad u square 2: 0.361384
Cycle 3:
   Number of active cells:   1280
   Number of degrees of freedom: 5185
 u square 3: 0.0549874  grad u square 3: 0.365557
Cycle 4:
   Number of active cells:   5120
   Number of degrees of freedom: 20609
 u square 4: 0.0553762  grad u square 4: 0.367069
Cycle 5:
   Number of active cells:   20480
   Number of degrees of freedom: 82177
 u square 5: 0.0555421  grad u square 5: 0.367712
Cycle 6:
   Number of active cells:   81920
   Number of degrees of freedom: 328193
 u square 6: 0.0556186  grad u square 6: 0.367986
Cycle 7:
   Number of active cells:   327680
   Number of degrees of freedom: 1311745
 u square 7: 0.0556576  grad u square 7: 0.368119


and the ratio 

[image: ask-8-1.JPG]
are:
u_square:1.84102.40381.31491.22871.1168
0.9720
grad_u_square:   2.06522.27661.46461.23361.28430.9891

And the results without adaptive mesh seem like approximate the exact 
results correctly. So I am confused why the results of "grad u square" with 
adaptive mesh are so strange? Maybe the reason is that we refinre mesh 
where the gradient of u is huge, so only "grad_u_square" can be very 
different with changing adaptive meshes? If so, how to deal with this 
problem? How to approximate   "grad u_ecaxt square" correctly with changing 
adaptive mesh? Maybe change the value smaller of scale of refine cells and 
coarsen cells, to make the change smaller?

Thank you very much!

Best,
Chucui
 

-- 
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] different results when compute integration with adaptive mesh

[chucui1016]

Hi all,

I have a very simple question but cannot think out:

I use step-6.cc to implement different adaptive mesh, and solve u with each 
mesh. Then I want to compute : 

[image: ask-7.JPG]
But there is a big gup on the results of the integration of the square of 
grad u:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
 u square 0: 0.0305621  grad u square 0: 0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
 u square 1: 0.0438936  grad u square 1: 1.31255
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
 u square 2: 0.0505902  grad u square 2: 3.584
Cycle 3:
   Number of active cells:   200
   Number of degrees of freedom: 921
 u square 3: 0.0536183  grad u square 3: 6.60495
Cycle 4:
   Number of active cells:   440
   Number of degrees of freedom: 2017
 u square 4: 0.0555872  grad u square 4: 19.7488
Cycle 5:
   Number of active cells:   956
   Number of degrees of freedom: 4425
 u square 5: 0.0562045  grad u square 5: 48.0924
Cycle 6:
   Number of active cells:   1916
   Number of degrees of freedom: 8993
 u square 6: 0.0572988  grad u square 6: 117.562
Cycle 7:
   Number of active cells:   3860
   Number of degrees of freedom: 18353
 u square 7: 0.0577537  grad u square 7: 258.867


And I don't know why the differences are so huge, and attach is my test 
code.

Thanks in advance!

Best,
Chucui



-- 
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) 2000 - 2016 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 at
 * the top level of the deal.II distribution.
 *
 * -

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


// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

#include 

#include 

#include 

// Finally, this is as in previous programs:
using namespace dealii;

template 
class Step6
{
public:
  Step6 ();
  ~Step6 ();

  void run ();

private:
  void setup_system ();
  void assemble_system ();
  void solve ();
  void refine_grid ();
  void output_results (const unsigned int cycle) const;

  Triangulation   triangulation;

  DoFHandler  dof_handler;
  FE_Qfe;

  ConstraintMatrix constraints;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;
  SparseMatrix volume_matrix;
  SparseMatrix gradient_matrix;

  Vector   solution;
  Vector   system_rhs;
};


template 
double coefficient (const Point &p)
{
  if (p.square() < 0.5*0.5)
return 20;
  else
return 1;
}

template 
Step6::Step6 ()
  :
  dof_handler (triangulation),
  fe (2)
{}

template 
Step6::~Step6 ()
{
  dof_handler.clear ();
}

template 
void Step6::setup_system ()
{
  dof_handler.distribute_dofs (fe);

  solution.reinit (dof_handler.n_dofs());
  system_rhs.reinit (dof_handler.n_dofs());


  constraints.clear ();
  DoFTools::make_hanging_node_constraints (dof_handler,
   constraints);


  VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction(),
constraints);


  constraints.close ();

  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler,
  dsp,
  constraints,
  /*keep_constrained_dofs = */ false);

  sparsity_pattern.copy_from(dsp);

  system_matrix.reinit (sparsity_pattern);
  volume_matrix.reinit (sparsity_pattern);
  gradient_matrix.reinit (sparsity_pattern);
}

template 
void Step6::assemble_system ()
{
  const QGauss  quadrature_formula(3);

  FEValues fe_values (fe, quadrature_formula,

[deal.II] Re: Implement of adaptive mesh refinement cannot work in time-dependent problem.

[chucui1016]

Dear Prof. Arndt,

Thank you very much!  
>when the solver fails to converge within 1000 iterations.
This is just because my solver code is :
  template 
  void StokesProblem::solve_preconditioner ()
  { 
SolverControl solver_control (1000, 1e-8*system_rhs.l2_norm());
SolverGMRES > gmres (solver_control);
gmres.solve (system_matrix, solution, system_rhs, 
PreconditionIdentity());//, preconditioner_H
  }

when I change "1000" above into larger number, like "10", the original 
code can run. And I rewrite the whole code, and the code can run now. 
Thankyou very much!

Best,
Chucui

在 2019年1月16日星期三 UTC+8下午7:42:28，Daniel Arndt写道：
>
> Chucui,
>
>
> Thanks for your quick answer! In my deal.ii (version 8.5.1), the code also 
>> fails in the 9th time step. Shall I need to reinstall my deal.ii into the 
>> newest version? 
>>
>  
> Normally, we only fix bugs in the library in the developer version which 
> is also very stable unless you are trying to use the newest features.
> However, I would expect that there is a way to write your code in such a 
> way that it also works with version 8.5 correctly.
>
>
>  
>
>>  For the shortest test code, please give me some time. Because I have 
>> deleted about 2000 lines. If I delete some other lines, the error 
>> information will change. 
>>
>  
> The faster you do that, the faster we can help you. It is totally up to 
> you. ;-) 
>
> 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: Implement of adaptive mesh refinement cannot work in time-dependent problem.

[chucui1016]

Dear Prof. Arndt,

Thanks for your quick answer! In my deal.ii (version 8.5.1), the code also 
fails in the 9th time step. Shall I need to reinstall my deal.ii into the 
newest version? 

 For the shortest test code, please give me some time. Because I have 
deleted about 2000 lines. If I delete some other lines, the error 
information will change. 

Thank you!

Best,
Chucui

在 2019年1月16日星期三 UTC+8下午6:28:27，Daniel Arndt写道：
>
> Chucui,
>
> with a recent developer version your code fails for me in the 9th time 
> step when the solver fails to converge within 1000 iterations.
> Hence, I can't reproduce the problem with a newer deal.II version.
>
> Please provide a minimal example by removing all functionality that is not 
> necessary to show the issue.
> It is simply not feasible to step through a code with almost 4000 lines of 
> code.
>
> If I had to guess, I would say that you are calling clear() on one of your 
> BlockSparseMatrix objects before the object is destroyed
> (without reinitializing the object in between) and this case is not 
> handled correctly in version 8.5.
>
> 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: Implement of adaptive mesh refinement cannot work in time-dependent problem.

[chucui1016]

Dear Prof. Wolfgang,

Thank you very much for your quick answer! I have thinked about your 
answers and debuged my code these days, and my question No.1 is still 
unsolved, and question No.2 has been solved:
 question No.1 : how to destroy sp efficiently? (unsolved)
As you say, I can destroy SparsityPattern by using the marked line below, 
which is " m.reinit(); "
>{ 
>   SparseMatrixm; 
>   SparsityPattern sp; 
>   ...fill sp, for example by copying from a DynamicSparsityPattern... 
>   m.reinit (sp);  // m now points to sp 
>   ...some more code... 

>   m.reinit();  // *** 
>}
But is I use the line directly, which means there isn't anything in the 
brackets:
{
BlockSparseMatrix system_matrix;
BlockSparsityPattern  sparsity_pattern_1;
[...]
sparsity_pattern_1.copy_from (dsp);
system_matrix.reinit (sparsity_pattern_1);
[...]
system_matrix.reinit ();
}

And the code cannot run, and the error is shown as :
error: no matching function for call to 
‘dealii::BlockSparseMatrix::reinit()’
 system_matrix.reinit ();

If I define another BlockSparsityPattern : sparsity_pattern_2; then 
{
BlockSparseMatrix system_matrix;
BlockSparsityPattern  sparsity_pattern_1;
BlockSparsityPattern  sparsity_pattern_2;
[...]
sparsity_pattern_1.copy_from (dsp);
sparsity_pattern_2.copy_from (dsp);
system_matrix.reinit (sparsity_pattern_1);
[...]
system_matrix.reinit (sparsity_pattern_2);
}
So the sparsity_pattern_1 is destroyed, but sparsity_pattern_2 cannot to be 
destroyed, there is a BlockSparsityPattern object cannot be destroyed in 
each time step, How to deal with it?


 question No.2 : what differences between triangulation and 
triangulation_active? (solved)
The definition of these two are:
GridGenerator::hyper_cube (triangulation_active, 0, 4.5);
triangulation_active.refine_global (8);
GridGenerator::flatten_triangulation (triangulation_active, 
triangulation);
By using GridGenerator::flatten_triangulation, the level of triangulation 
is 0, i.e. the coarsest mesh, which is used when define a periodic boundary 
conditions. So triangulation.n_levels() is always 0. and cannot use 
triangulation when determine which cell to be coarsen or refined:
if (triangulation.n_levels() > max_grid_level) 
>for (typename Triangulation::active_cell_iterator 
> cell = triangulation.begin_active(max_grid_level); 
> cell != triangulation.end(); ++cell) 
>  cell->clear_refine_flag (); 

So if we not only want to use triangulation.n_levels() like above,but also 
want to use the coarsest mesh to define periodic boundary conditions,  we 
should learn from step-45, like this:
GridGenerator::hyper_cube (triangulation, 0, 4.5);


for (typename Triangulation::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f::faces_per_cell; 
++f)
{
   if (cell->face(f)->at_boundary())
   {  
  if (std::fabs(cell->face(f)->center()(0) - (4.5)) 
< 1e-12)
  
  cell->face(f)->set_boundary_id (1);
  else if(std::fabs(cell->face(f)->center()(1) - 
(0)) < 1e-12)
   
   cell->face(f)->set_boundary_id (2);
   else if(std::fabs(cell->face(f)->center()(1) 
- (4.5)) < 1e-12)
 
 cell->face(f)->set_boundary_id (3);
 else 
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
  {
  cell->face(f)->set_boundary_id 
(0);
  }  
   } 
   }
}


triangulation.refine_global (3);

Hope the answer of question No.1. Thanks in advance. 

Best,
Chucui




>

-- 
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: Implement of adaptive mesh refinement cannot work in time-dependent problem.

[chucui1016]

 Dear Prof. Wolfgang,

Thank you very much for your quick answer! And I have tried what you say, 
and I have some questions still:
 
1.

> In your case, a 
> SparsityPattern is still used by some other object but is being destroyed. 
> You'll have to find out what other object that is -- likely a SparseMatrix 
> object. You'll first have to break that dependence before you destroy the 
> sparsity pattern. 
>
 I had read : 
https://groups.google.com/forum/#!searchin/dealii/Object$20of$20class$20N6dealii15SparsityPatternE$20is$20still$20used$20by$201$20other$20objects%7Csort:date/dealii/8Sqkbo1n5Jg/O6R2Vpya4hwJ
And I think what you means is the same as it above beforehand, but maybe 
not now. The dsp (BlockDynamicSparsityPattern) will be destroyed and cannot 
exist in all lifetime of  BlockSparseMatrix, so we define 
: BlockSparsityPattern  sparsity_pattern at the top, (this is what 
step-22 does) and let it copy dsp
But now, we know that BlockSparsityPattern  sparsity_pattern also will 
be destroyed and cannot exist in all lifetime of  BlockSparseMatrix.  So I 
don't know how to correct it.
And I cannot implement what you said "You'll first have to break that 
dependence before you destroy the sparsity pattern".  

2. As your suggestions, I change my code :
  template 
  void
  StokesProblem::refine_mesh (const unsigned int max_grid_level, 
   const unsigned int min_grid_level)
  {
Vector estimated_error_per_cell 
(triangulation_active.n_active_cells());
FEValuesExtractors::Scalar phase(0);   
KellyErrorEstimator::estimate (dof_handler,
QGauss(degree+1),
typename FunctionMap::type(),
old_solution,
estimated_error_per_cell,
fe.component_mask(phase)); 
GridRefinement::refine_and_coarsen_fixed_number (triangulation_active,

 estimated_error_per_cell,
 0.2,0.1);  
std::cout << "triangulation.n_levels" << triangulation_active.n_levels() << 
std::endl;  
   
   if (triangulation_active.n_levels() > max_grid_level)
  for (typename Triangulation::active_cell_iterator
   cell = triangulation.begin_active(max_grid_level);
   cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();   
std::cout << "clear-refine" << triangulation_active.n_levels() << 
std::endl;
   if (triangulation_active.n_levels() < min_grid_level)
  for (typename Triangulation::active_cell_iterator
   cell = triangulation.end_active(min_grid_level);
   cell != triangulation.end(); ++cell)
cell->clear_coarsen_flag ();   
std::cout << "clear-coarsen" << triangulation_active.n_levels() << 
std::endl;  
triangulation_active.execute_coarsening_and_refinement ();
std::cout << "c-and-r" << triangulation_active.n_levels() << std::endl; 
  }

Maybe I need to chage all "triangulation" in my originial 
StokesProblem::refine_mesh 
(const unsigned int max_grid_level, const unsigned int min_grid_level) into 
"triangulation_active", this function can run(but i don't know if it is 
right), and for implement periodic boundary condition, I use the "
triangulation" and  "triangulation_active" :
for (unsigned int cycle=0; cycle::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f::faces_per_cell; 
++f)
{
   if (cell->face(f)->at_boundary())
   {  
  if (std::fabs(cell->face(f)->center()(0) - (4.5)) 
< 1e-12)
  
  cell->face(f)->set_boundary_id (1);
  else if(std::fabs(cell->face(f)->center()(1) - 
(0)) < 1e-12)
   
   cell->face(f)->set_boundary_id (2);
   else if(std::fabs(cell->face(f)->center()(1) 
- (4.5)) < 1e-12)
 
 cell->face(f)->set_boundary_id (3);
 else 
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
  {
  cell->face(f)->set_boundary_id 
(0);
  }  
   } 
   }
}
 
 }
I don't know what I should use in whole code is  "triangulation"

[deal.II] Implement of adaptive mesh refinement cannot work in time-dependent problem.

[chucui1016]

Dear All,

I have 2 questions when implementing of adaptive mesh refinement in 
time-dependent problem, and the attached is my test code: test_for_ask.cc

1. I use the form to update my mesh:
[...]
timestep_number = 1;
setup_dofs ();
[...]
for (timestep_number=2; timestep_number<=100; ++timestep_number) {
[...]
refine_mesh (max_grid_level, min_grid_level);
setup_dofs ();
[...]
}

But starting from timestep_number=2, the setup_dof () cannot work. the 
error is:

An error occurred in line <128> of file 
 in function
void dealii::Subscriptor::check_no_subscribers() const
The violated condition was: 
counter == 0
Additional information: 
Object of class N6dealii15SparsityPatternE is still used by 1 other 
objects.

(Additional information: 
  from Subscriber SparseMatrix)

See the entry in the Frequently Asked Questions of deal.II (linked to from 
http://www.dealii.org/) for a lot more information on what this error means 
and how to fix programs in which it happens.

Stacktrace:
---
#0  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::Subscriptor::check_no_subscribers() const
#1  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::Subscriptor::~Subscriptor()
#2  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::SparsityPattern::~SparsityPattern()
#3  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::SparsityPattern::~SparsityPattern()
#4  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::BlockSparsityPatternBase::reinit(unsigned 
int, unsigned int)
#5  //home/y/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: 
dealii::BlockSparsityPattern::copy_from(dealii::BlockDynamicSparsityPattern 
const&)
#6  ./Problem13-2: Step22::StokesProblem<2>::setup_dofs()
#7  ./Problem13-2: Step22::StokesProblem<2>::run()
#8  ./Problem13-2: main


make[3]: *** [CMakeFiles/run] Aborted (core dumped)
make[2]: *** [CMakeFiles/run.dir/all] Error 2
make[1]: *** [CMakeFiles/run.dir/rule] Error 2
make: *** [run] Error 2


And my setup_dofs () is written by learning from step-22, which means it 
use BlockSparsityPattern:
[...]
system_matrix.clear ();
system_matrix_preconditioner.clear ();
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);   
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
   

dsp.collect_sizes();

DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
}
system_matrix.reinit (sparsity_pattern);
system_matrix_preconditioner.reinit (sparsity_pattern);
[...]

 I had tried several methods all cannot work:(if I set another 
BlockSparsityPattern sparsity_pattern_new, and use it in timestep_number = 
2, it can go, but stop at timestep_number = 3; or I use 
sparsity_pattern_pointer
= std_cxx11::shared_ptr::type>(new 
typename BSP_solution::type()); it also doesn't work, where:
  template 
  struct BSP_solution;
  template <>
  struct BSP_solution<2>
  {
typedef BlockSparsityPattern type;
  };



 and now I don't know how to correct it. Can anyone give me a help?


2. My refine_mesh() is: 
  template 
  void
  StokesProblem::refine_mesh (const unsigned int max_grid_level, 
   const unsigned int min_grid_level)
  {
Vector estimated_error_per_cell (triangulation.n_active_cells());
FEValuesExtractors::Scalar phase(0);   
KellyErrorEstimator::estimate (dof_handler,
QGauss(degree+1),
typename FunctionMap::type(),
old_solution,
estimated_error_per_cell,
fe.component_mask(phase)); 
GridRefinement::refine_and_coarsen_fixed_number (triangulation,

 estimated_error_per_cell,
 0.2,0.1);  

   if (triangulation.n_levels() > max_grid_level)
  for (typename Triangulation::active_cell_iterator
   cell = triangulation.begin_active(max_grid_level);
   cell != triangulation.end(); ++cell)
cell->clear_refine_flag (); 
/*  this part for clear_coarsen_flag() cannot work */  
  
   if (triangulation.n_levels() < min_grid_level)
  for (typename Triangulation::active_cell_iterator
   cell = triangulation.begin_active(min_grid_level);
   cell != triangulation.end(); ++cell)
cell->clear_coarsen_flag ();
 
triangulation.

Re: [deal.II] How to compute the cellwise error of the finite element solution without exact solution?

[chucui1016]

Dear Prof. Bangerth,

Thank you very much! 
FEFieldFunction fe_function_finest (dof_handler_finest, solution_finest
);
And I can use the fe_function_finest in function:  VectorTools::
integrate_difference 

And now I have another question: I want to compute numerical solutions in 
different widths of meshgrid, so I need different triangulations, 
dof_handlers, solution vectors. And I need to compute the finest one at 
first because I need to use the solution_finest as the 'exact solution'. 
The question is how to coarse the meshgrid after getting the 
solution_finest? If I want to refine the mesh, I will use 
'triangulation.refine(1)'. But I don't know how to coarse the mesh. Or 
maybe I have to define different triangulations, dof_handlers at first, and 
compute different solution vectors. If so, the code will be very long. 
Thank you very much!

Best,
Chucui

>
> You can't do this directly -- the function just wants an object that can 
> be evaluated at arbitrary points. Passing a vector by itself will not 
> help: a vector is just a bunch of numbers, but it does not by itself 
> explain how it is to be interepreted. Indeed, a vector by itself is 
> meaningless unless you also specify a DoFHandler that makes the 
> connection between nodal values, where these nodes are located in 
> physical space, and how the shape functions are defined. 
>
> But there is a class that combines all of this knowledge and that then 
> represents a finite element field (=DoFHandler + Vector) as a function: 
> FEFieldFunction. Take a look at that class, as it implements the 
> interface of a Function object that can be passed in the place you point 
> out. It is not particularly fast, but should allow you to do what you 
> want. 
>
> Of course, you have to make sure that both the 'dof' and 'fe_function' 
> arguments to the VectorTools::integrate_difference() function match 
> (i.e., you can't just save a solution vector, and then refine the 
> triangulation and thereby change the DoFHandler), as well as the 
> DoFHandler and Vector objects you pass to FEFieldFunction. 
>
> 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] Re: How to restart a time dependent code from its breakpoint?

[chucui1016]

Dear Jean-Paul and dear Stephen,

Thank you for your reply! As the examples you set, and an answer I 
found:  See "Use for serialization" in 
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_developer_doxygen_deal.II_classparallel-5F1-5F1distributed-5F1-5F1SolutionTransfer.html&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=88KR9zh1f_gaVgWHpZwwXXlaOCduKgZ1yzsuERztdmU&s=Q549G3WmW1DtS06Yq9GyH98gIg8Rl44Ezac7I1Fc1Tc&e=
 

 
these functions such as  ' Triangulation::save()', 'Triangulation::load()', 
void SolutionTransfer< dim, VectorType, DoFHandlerType 
>::prepare_serialization ( const VectorType &  in ) 
 void SolutionTransfer< dim, VectorType, DoFHandlerType >::deserialize ( 
VectorType 
&  in )
are all used in the parallel examples.
But my code is not use any class and headfile for parallel one. Shall I 
rewrite my code into a parallel one? Thank you very much!

Best,
Chucui 


-- 
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] How to restart a time dependent code from its breakpoint?

[chucui1016]

Dear Jean-Paul,

Thank you for your reply! As the examples you set, and an answer I 
found:  See "Use for serialization" in 
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_developer_doxygen_deal.II_classparallel-5F1-5F1distributed-5F1-5F1SolutionTransfer.html&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=88KR9zh1f_gaVgWHpZwwXXlaOCduKgZ1yzsuERztdmU&s=Q549G3WmW1DtS06Yq9GyH98gIg8Rl44Ezac7I1Fc1Tc&e=
 

 
these functions such as  ' Triangulation::save()', 'Triangulation::load()', 
void SolutionTransfer< dim, VectorType, DoFHandlerType 
>::prepare_serialization ( const VectorType & in ) 
 void SolutionTransfer< dim, VectorType, DoFHandlerType >::deserialize ( 
VectorType 
& in )
are all used in the parallel examples.
But my code is not use any class and headfile for parallel one. Shall I 
rewrite my code into a parallel one? Thank you very much!

Best,
Chucui   

在 2018年11月5日星期一 UTC+8下午5:01:06，Jean-Paul Pelteret写道：
>
> Dear Chucui,
>
> I don’t use this functionality myself, so I can only point you to the 
> test suite 
>  where 
> its functionality is tested. Here 
>  is 
> one example where a distributed triangulation is saved and loaded, and 
> here 
>  is 
> one where a distributed triangulation + solution is saved and loaded with a 
> different number of MPI processes.
>
> Best,
> 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] How to compute the cellwise error of the finite element solution without exact solution?

[chucui1016]

Dear all,

If I don't have the exact solution, but I want to compute the L2 form of 
u_{h}-u_{finest} (i.e. see u_{finest} as the exact solution, and u_{finest} 
is the numerical solution with the finest mesh grid among the series we 
want to test, for example, the width of the unit meshgrid is h = 1/2^k, 
k=6,7,8,9,10. So the u_{finest}=u_{h=1/2^10}). Does deal.ii have a function 
like
template

void VectorTools::integrate_difference ( const Mapping 
< dim, 
spacedim > & mapping,
const DoFHandler 
< dim
, spacedim > & dof,
const InVector & fe_function,
const Function 
< 
spacedim, double > & exact_solution,
OutVector & difference,
const Quadrature 
< 
dim > & q,
const NormType 

 
& norm,
const Function 
< 
spacedim, double > * weight = nullptr,
const double exponent = 2. 
)


and doesn't use 

const Function 
<
 spacedim, double > & exact_solution,
but use

const InVector & fe_function,
at this location of input parameters. Thanks in advance!

Best, 
Chucui

-- 
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] How to restart a time dependent code from its breakpoint?

[chucui1016]

Dear Jean-Paul,

Thanks for your reply! And I am interested in 'the serialization 
capabilities', can you give me a little example or code to show me how it 
work. For me the ASPECT library 

 is 
an entile new one, I cannot find the example that shows the serialization 
capabilities exactly, so I need some help. And is the point of  'the 
serialization capabilities' a fundamental point of C++? Does it only work 
at some projects such as  the ASPECT library 

 you 
mentioned above? Thank you very much!

Best,
Chucui

在 2018年11月5日星期一 UTC+8下午3:36:59，Jean-Paul Pelteret写道：
>
> Dear Chucui,
>
> The Triangulation 
> 
> , DoFHandler 
> 
>  and Vector 
> 
>  classes 
> support serialization, that is to say that they have built-in save() / 
> load() functions. The GridIn::read_vtk() 
> 
>  function 
> can supposedly construct a triangulation from a solution vector, but I’m 
> not sure if this is working at the moment (see this GitHub issue 
> ). I don’t think that there 
> exists a mechanism to read in the numerical data from VTK file as of yet. 
> So I would suggest that you utilise the serialization capabilities, which 
> should allow a complete restart of a simulation from any tilmestep. I 
> understand that this is what the ASPECT library 
>  does.
>
> Although it does not directly help you restart a simulation using VTK 
> data, I hope that you find this information useful.
>
> Best,
> 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: How to restart a time dependent code from its breakpoint?

[chucui1016]

Dear Jonathan,

Thank you very much! Firstly, My vtk file's names are allowed yo have more 
than three digits. My code of writting vtk file is:
  template 
  void MyProblem::output_results_U (const unsigned int 
timestep_number)  const
  { 
DataOut data_out;

data_out.attach_dof_handler (dof_handler_U);
data_out.add_data_vector (old_U, "U");

data_out.build_patches (degree+3);

std::ostringstream filename;
filename << "U-"
 << Utilities::int_to_string(timestep_number)
 << ".vtk";

std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
  }
And I have got VTK files such as 'U-6400.vtk'

I'm sorry for my discription without codes. For example, the system I want 
to solve is just like:

[image: ol.JPG]

Then the weak form in each time step is:

[image: 哦了.JPG]

The code to assemble system matrix and right hand side is:
  template 
  void MyProblem::assemble_system_U (const Vector old_U)
  {
system_matrix_U=0;
system_rhs_U = 0;
QGauss   quadrature_formula(degree+2);

FEValues fe_U_values (fe_U, quadrature_formula,
   update_values| update_gradients |
   update_quadrature_points | update_hessians | 
update_JxW_values); 
   
const unsigned int   dofs_per_cell_U   = fe_U.dofs_per_cell;
const unsigned int   n_q_points  = quadrature_formula.size();
Vector   local_rhs (dofs_per_cell_U);
FullMatrix   local_matrix (dofs_per_cell_U, dofs_per_cell_U);
std::vector local_dof_indices 
(dofs_per_cell_U);
   
std::vector old_U_values(n_q_points);
  
typename DoFHandler::active_cell_iterator
cell_U = dof_handler_U.begin_active(),
endc_U = dof_handler_U.end();
for (; cell_U!=endc_U; ++cell_U)
  {
local_rhs = 0;
local_matrix=0;
fe_U_values.reinit (cell_U);


fe_U_values.get_function_values (old_U, old_U_values);

for (unsigned int q=0; qget_dof_indices (local_dof_indices);
constraints_U.distribute_local_to_global(local_matrix, local_rhs, 
local_dof_indices, system_matrix_U, system_rhs_U);
  }  
  }
Then solve new_U i.e. U_{n+1}, and update the U vector.

template 
  void MyProblem::solve_U ()
  { 

SparseDirectUMFPACK  A_direct;
A_direct.initialize(system_matrix_U);
A_direct.vmult (new_U, system_rhs_U);
  } 
  template 
  void MyProblem::run ()
  {
[...]
VectorTools::project (dof_handler_U, constraints_U, 
QGauss(degree+2),
InitialValues_U(),
old_U);
int timestep_number = 1;
 for (; time<=1.0; time+=time_step, ++timestep_number)  
 {
assemble_system_U (old_U);   
solve_U (); 
constraints_U.distribute (new_U); 
 old_U = new_U; 
output_results_U (timestep_number);}}

But if the code is broken suddenly ,at time step n=1000 for example, how to 
restart it from time step t_n = 1001, not t_n=1? And in each time step, I 
have a output VTK file of solution U (U-n.vtk), but in the code I need 
‘const Vector old_U’ in fact. How to convert a U-1000.vtk file into 
a 'const Vector old_U', then use it for 'assemble_system_U (const 
Vector old_U)'?  Thank you very much!

Best,
Chucui 

在 2018年11月5日星期一 UTC+8下午1:25:58，mrjonm...@gmail.com写道：
>
> If you're using the timestep number in the file names of your vtk files, 
> are you formatting the filename to allow more than three digits? (e.g. 
> solution-999.vtk versus solution-0999.vtk.)
>
> It's hard to tell without seeing any of your code.
> Cheers,
> Jonathan 
>
>

-- 
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] How to restart a time dependent code from its breakpoint?

[chucui1016]

Hi All,

I am running a code where the system I want to solve is time-dependent, 
which means:

[image: 快快快.JPG]

So, I want to solve u_{n+1} by using u_{n} in each time step n.
If my code is broken suddenly ,at time step n=1000 for example, how to 
restart it from time step t_n = 1001, not t_n=1? And in each time step, I 
have a output VTK file of solution u.

Thank you in advance!

Best,
Chucui

-- 
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: After "make install", there are some Errors when install Deal.II-9.0.1

[chucui1016]

Dear Wolfgang,

Thanks for your advice, and Matthias's. Thank you very much!

>
>
> I don't think any of us has seen this error before. Can you try what 
> happens if you delete the build directory, re-create it, and then just 
> use the minimal command line 
>cmake /vol7/... 
> ? Make things as simple as possible if you don't know what causes a 
> particular problem. If the problem goes away, then repeat by adding one 
> flag at a time until the problem re-appears -- at which time you know 
> what change must have caused it. 
>
 I download a new dealii-9.0.1.tar.gz and build a new file folder:
/vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/ 
source directory:/vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/src
and the command line :(in the directory: 
/vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/src/build)
$ cmake /vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/src
$ make -j2
and it shows:
In file included from 
/vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/src/source/base/flow_function.cc:18:0:
/vol7/home/zhaoyucan1016/reinstall2-dealii-9.0.1/src/include/deal.II/base/flow_function.h:59:13:
 
internal compiler error: in use_thunk, at cp/method.c:338
 virtual ~FlowFunction() = default;
 ^
0x595884 use_thunk(tree_node*, bool)
../.././gcc/cp/method.c:338
0x5a078c emit_associated_thunks(tree_node*)
../.././gcc/cp/semantics.c:3763
0x5a0a47 expand_or_defer_fn(tree_node*)
../.././gcc/cp/semantics.c:3890
0x5b83a5 maybe_clone_body(tree_node*)
../.././gcc/cp/optimize.c:428
0x5a083d expand_or_defer_fn_1(tree_node*)
../.././gcc/cp/semantics.c:3814
0x5a0a18 expand_or_defer_fn(tree_node*)
../.././gcc/cp/semantics.c:3884
0x5969ef synthesize_method(tree_node*)
../.././gcc/cp/method.c:809
0x5434aa mark_used(tree_node*)
../.././gcc/cp/decl2.c:4677
0x4f01cc build_over_call
../.././gcc/cp/call.c:7055
0x4ed678 build_new_method_call_1
../.././gcc/cp/call.c:7715
0x4ed678 build_new_method_call(tree_node*, tree_node*, vec**, tree_node*, int, tree_node**, int)
../.././gcc/cp/call.c:7785
0x4edf1e build_special_member_call(tree_node*, tree_node*, vec**, tree_node*, int, int)
../.././gcc/cp/call.c:7352
0x58e0ff expand_cleanup_for_base
../.././gcc/cp/init.c:1217
0x5925f1 expand_cleanup_for_base
../.././gcc/cp/init.c:1112
0x5925f1 emit_mem_initializers(tree_node*)
../.././gcc/cp/init.c:1097
0x56a47c cp_parser_mem_initializer_list
../.././gcc/cp/parser.c:11731
0x56a47c cp_parser_ctor_initializer_opt
../.././gcc/cp/parser.c:11642
0x56a47c cp_parser_ctor_initializer_opt_and_function_body
../.././gcc/cp/parser.c:17848
0x56aacf cp_parser_function_definition_after_declarator
../.././gcc/cp/parser.c:21843
0x56b6e3 cp_parser_function_definition_from_specifiers_and_declarator
../.././gcc/cp/parser.c:21764
Please submit a full bug report,
with preprocessed source if appropriate.
Please include the complete backtrace with any bug report.
See  for instructions.
make[2]: *** [source/base/CMakeFiles/obj_base_debug.dir/flow_function.cc.o] 
Error 1
make[1]: *** [source/base/CMakeFiles/obj_base_debug.dir/all] Error 2
make[1]: *** Waiting for unfinished jobs
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_ilu.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_matrix_ez.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_mic.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_vanka.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparsity_pattern.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparsity_tools.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/swappable_vector.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/tridiagonal_matrix.cc.o
[ 42%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/vector.cc.o
[ 43%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/vector_memory.cc.o
[ 43%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/vector_view.cc.o
[ 43%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_matrix.cc.o
[ 43%] Building CXX object 
source/lac/CMakeFiles/obj_lac_debug.dir/sparse_matrix_inst2.cc.o
[ 43%] Built target obj_lac_debug
make: *** [all] Error 2

and  I download a new dealii-9.0.1.tar.gz and build a new file 
folder: /vol7/home/zhaoyucan1016/reinstall-dealii-9.0.1
source directory:/vol7/home/zhaoyucan1016/reinstall-dealii-9.0.1/src
and the command line :(in the directory: 
/vol7/home/zhaoyucan1016/reinstall-dealii-9.0.1/src/build)
$ cmake /vol7/home/zhaoyucan1016/reinstall-dealii-9.0.1/src 
-DCMAKE_INSTALL_PREFIX=/vol7/home/zhaoyucan1016/reinstall-dealii-9.0.1/deal.II 
-DDEAL_II_WITH_MPI=ON -DMPI_DIR=/usr/local/mpi3 
-DCMAKE_CXX_COMPILER=/usr/local/mpi3/bin/mpicxx 
-DCMAKE_Fortran_COMPILER=/usr/local/mpi3/bin/mpifort 
-DDEAL_II_FORCE_BUNDLED_BOOST=ON -DDEAL_II_ALLOW_AUTODETECTION=ON
$make -j2 deal_II
$make 

[deal.II] Re: After "make install", there are some Errors when install Deal.II-9.0.1

[chucui1016]

Dear Matthias,

The newest result are still the same. After "make install" it shows:

[ 82%] Built target obj_rol_debug
Scanning dependencies of target step-3.release
[ 82%] Building CXX object 
examples/CMakeFiles/step-3.release.dir/step-3/step-3.cc.o
[ 82%] Linking CXX executable ../bin/step-3.release
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 1>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<3, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<3, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 1>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 2>::StraightBoundary()'
make[2]: *** [bin/step-3.release] Error 1
make[1]: *** [examples/CMakeFiles/step-3.release.dir/all] Error 2
make: *** [all] Error 2

And  I am on the 
directory: /vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1

My command lines are
$cmake /vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1  
-DDEAL_II_WITH_MPI=ON  -DCMAKE_C_COMPILER=/usr/local/mpi3/bin/mpicc  
-DCMAKE_CXX_COMPILER=/usr/local/mpi3/bin/mpicxx  
-DCMAKE_Fortran_COMPILER=/usr/local/mpi3/bin/mpifort 
-DCMAKE_INSTALL_PREFIX=/vol7/home/zhaoyucan1016/dealii-9.0.1/build 
-DMPI_DIR=/usr/local/mpi3 
-DP4EST_DIR=/vol7/home/zhaoyucan1016/p4est/p4est-2.0 
-DPETSC_DIR=/vol7/home/zhaoyucan1016/petsc-3.7.7/petsc-3.7.7 
-DHDF5_DIR=/usr/lib64/ 
-DMETIS_DIR=/vol7/home/zhaoyucan1016/metis-5.1.0/metis-5.1.0 
-DDEAL_II_WITH_METIS=ON -DDEAL_II_WITH_HDF5=OFF -DDEAL_II_WITH_PETSC=ON 
-DDEAL_II_WITH_P4EST=OFF -DDEAL_II_WITH_MPI=ON 
-DDEAL_II_FORCE_BUNDLED_BOOST=ON-DDEAL_II_ALLOW_AUTODETECTION=ON 
$make -j2 deal.II
$make install

And under the directory: /vol7/home/zhaoyucan1016/dealii-9.0.1/build
There are only : CMakeCache.txt  CMakeFiles (Folder)

Hope these informations can help.

Best,
Chucui


-- 
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 configuration:
#CMAKE_BUILD_TYPE:   DebugRelease
#BUILD_SHARED_LIBS:  ON
#CMAKE_INSTALL_PREFIX:   /vol7/home/zhaoyucan1016/dealii-9.0.1/build
#CMAKE_SOURCE_DIR:   
/vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1
#(version 9.0.1)
#CMAKE_BINARY_DIR:   
/vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1
#CMAKE_CXX_COMPILER: Intel 15.0.1.20141023 on platform Linux x86_64
#/usr/local/mpi3/bin/mpicxx
#CMAKE_C_COMPILER:   /usr/local/mpi3/bin/mpicc
#CMAKE_Fortran_COMPILER: /usr/local/mpi3/bin/mpifort
#CMAKE_GENERATOR:Unix Makefiles
#
#  Base configuration (prior to feature configuration):
#DEAL_II_CXX_FLAGS:-fpic -ansi -w2 -diag-disable=remark 
-wd21 -wd68 -wd135 -wd175 -wd177 -wd191 -wd193 -wd279 -wd327 -wd383 -wd981 
-wd1418 -wd1478 -wd1572 -wd2259 -wd2536 -wd3415 -wd15531 -wd111 -wd128 -wd185 
-wd186 -wd280 -qopenmp-simd -std=c++11
#DEAL_II_CXX_FLAGS_RELEASE:-O2 -no-ansi-alias -ip -funroll-loops 
-no-vec
#DEAL_II_CXX_FLAGS_DEBUG:  -O0 -g -gdwarf-2 -grecord-gcc-switches
#DEAL_II_LINKER_FLAGS: -Wl,--as-needed -shared-intel -qopenmp 
-rdynamic
#DEAL_II_LINKER_FLAGS_RELEASE: 
#DEAL_II_LINKER_FLAGS_DEBUG:   
#DEAL_II_DEFINITIONS:  
#DEAL_II_DEFINITIONS_RELEASE:  
#DEAL_II_DEFINITIONS_DEBUG:DEBUG
#DEAL_II_USER_DEFINITIONS: 
#DEAL_II_USER_DEFINITIONS_REL: 
#DEAL_II_USER_DEFINITIONS_DEB: DEBUG
#DEAL_II_INCLUDE_DIRS  
#DEAL_II_USER_INCLUDE_DIRS:
#DEAL_II_BUNDLED_INCLUDE_DIRS: 
#DEAL_II_LIBRARIES:m
#DEAL_II_LIBRARIES_RELEASE:
#DEAL_II_LIBRARIES_DEBUG:  
#DEAL_II_COMPILER_VECTORIZATION_LEVEL: 1
#
#  Configured Features (

Re: [deal.II] Re: After "make install", there are some Errors when install Deal.II-9.0.1

[chucui1016]

Dear Matthias,

I am sorry to delete the originial file and configure it again, but 
everything seems not change. And attached is the newest detailed.log. Thank 
you very much!

Best,
Chucui

在 2018年10月19日星期五 UTC+8上午11:53:35，Matthias Maier写道：
>
>
> On Thu, Oct 18, 2018, at 17:48 CDT, David Wells  > wrote: 
>
> > Hi Chucui, 
> > 
> > StraightBoundary was removed in deal.II 9.0; it should not show up in 
> any 
> > object files so these errors don't make any sense. Did you use a clean 
> > build directory? 
>
> Not quite. It got completely removed in the current development version. 
>
> In version 9.0.1 it is present in the library. What is surprising though 
> is the fact that we completely instantiate the class: 
>
> grid/tria_boundary.inst.in (on branch dealii-9.0 
>  18 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : 
>  SPACE_DIMENSIONS) 
>  19 { 
> 
>  20 #if deal_II_dimension <= deal_II_space_dimension   
>   
>  21 template class Boundary deal_II_space_dimension>; 
>  22 template class StraightBoundary deal_II_space_dimension>; 
>  23 #endif 
>   
>  24 } 
>
> yet, Chucui observes an undefined reference to the constructor. Which is 
> a bit puzzling. 
>
> Chucui, can you please post your full `detailed.log` please? Further, 
> can you please try to recompile again by (a) unpacking dealii-9.0.1 to a 
> new directory, and configure and compile in a new build directory? 
>
> 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.
For more options, visit https://groups.google.com/d/optout.
###
#
#  deal.II configuration:
#CMAKE_BUILD_TYPE:   DebugRelease
#BUILD_SHARED_LIBS:  ON
#CMAKE_INSTALL_PREFIX:   /vol7/home/zhaoyucan1016/dealii-9.0.1/build
#CMAKE_SOURCE_DIR:   
/vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1
#(version 9.0.1)
#CMAKE_BINARY_DIR:   
/vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1
#CMAKE_CXX_COMPILER: Intel 15.0.1.20141023 on platform Linux x86_64
#/usr/local/mpi3/bin/mpicxx
#CMAKE_C_COMPILER:   /usr/local/mpi3/bin/mpicc
#CMAKE_Fortran_COMPILER: /usr/local/mpi3/bin/mpifort
#CMAKE_GENERATOR:Unix Makefiles
#
#  Base configuration (prior to feature configuration):
#DEAL_II_CXX_FLAGS:-fpic -ansi -w2 -diag-disable=remark 
-wd21 -wd68 -wd135 -wd175 -wd177 -wd191 -wd193 -wd279 -wd327 -wd383 -wd981 
-wd1418 -wd1478 -wd1572 -wd2259 -wd2536 -wd3415 -wd15531 -wd111 -wd128 -wd185 
-wd186 -wd280 -qopenmp-simd -std=c++11
#DEAL_II_CXX_FLAGS_RELEASE:-O2 -no-ansi-alias -ip -funroll-loops 
-no-vec
#DEAL_II_CXX_FLAGS_DEBUG:  -O0 -g -gdwarf-2 -grecord-gcc-switches
#DEAL_II_LINKER_FLAGS: -Wl,--as-needed -shared-intel -qopenmp 
-rdynamic
#DEAL_II_LINKER_FLAGS_RELEASE: 
#DEAL_II_LINKER_FLAGS_DEBUG:   
#DEAL_II_DEFINITIONS:  
#DEAL_II_DEFINITIONS_RELEASE:  
#DEAL_II_DEFINITIONS_DEBUG:DEBUG
#DEAL_II_USER_DEFINITIONS: 
#DEAL_II_USER_DEFINITIONS_REL: 
#DEAL_II_USER_DEFINITIONS_DEB: DEBUG
#DEAL_II_INCLUDE_DIRS  
#DEAL_II_USER_INCLUDE_DIRS:
#DEAL_II_BUNDLED_INCLUDE_DIRS: 
#DEAL_II_LIBRARIES:m
#DEAL_II_LIBRARIES_RELEASE:
#DEAL_II_LIBRARIES_DEBUG:  
#DEAL_II_COMPILER_VECTORIZATION_LEVEL: 1
#
#  Configured Features (DEAL_II_ALLOW_BUNDLED = ON, DEAL_II_ALLOW_AUTODETECTION 
= ON):
#  ( DEAL_II_WITH_64BIT_INDICES = OFF )
#  ( DEAL_II_WITH_ADOLC = OFF )
#DEAL_II_WITH_ARPACK set up with external dependencies
#ARPACK_LINKER_FLAGS = 
#ARPACK_LIBRARIES = 
/usr/lib64/libarpack.so;/usr/lib64/libopenblas.so;mpifort;mpi;rt;glex;ifcore;imf;m;ipgo;irc;pthread;svml;gcc;gcc_s;irc_s;dl;c
#  ( DEAL_II_WITH_ASSIMP = OFF )
#DEAL_II_WITH_BOOST set up with bundled packages
#BOOST_BUNDLED_INCLUDE_DIRS = 
/vol7/home/zhaoyucan1016/dealii-9.0.1/dealii-9.0.1/bundled/boost-1.62.0/include
#BOOST_LIBRARIES = rt
#  ( DEAL_II_WITH_CUDA = OFF )
#  ( DEAL_II_WITH_CXX14 = OFF )
#  ( DEAL_II_WITH_CXX17 = OFF )
#  ( DEAL_II_WITH_GMSH = OFF )
#DEAL_II_WITH_GSL set up with external dependencies
#GSL_VERSION = 1.13
#GSL_INCLUDE_DIRS = /usr/include
#GSL_USER_INCLUDE_DIRS = /usr/include
# 

[deal.II] Re: After "make install", there are some Errors when install Deal.II-9.0.1

[chucui1016]

Dear David，

I am under the directory: ../dealii-9.0.1/dealii-9.0.1
and use the command line  in "cmake"

$cmake ../dealii-9.0.1/dealii-9.0.1 (which is the file that was unziped) 
-DCMAKE_INSTALL_PREFIX=../dealii-9.0.1/build  (which is the new and empty)

and this has errors;

Best,
Chucui
 

在 2018年10月19日星期五 UTC+8上午6:48:58，David Wells写道：
>
> Hi Chucui,
>
> StraightBoundary was removed in deal.II 9.0; it should not show up in any 
> object files so these errors don't make any sense. Did you use a clean 
> build directory?
>
> Thanks,
> David Wells
>
>
>>

-- 
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] After "make install", there are some Errors when install Deal.II-9.0.1

[chucui1016]

Dear all, 

I want to reinstall my deal.ii because I need to run some codes in parallel 
just like Step-40, step-55 and Step-32, and I choose deal.ii-9.0.1 
The command lines I use are: 

$cmake ../dealii-9.0.1/dealii-9.0.1  -DDEAL_II_WITH_MPI=ON  
-DCMAKE_C_COMPILER=/usr/local/mpi3/bin/mpicc  
-DCMAKE_CXX_COMPILER=/usr/local/mpi3/bin/mpicxx  
-DCMAKE_Fortran_COMPILER=/usr/local/mpi3/bin/mpifort 
-DCMAKE_INSTALL_PREFIX=../dealii-9.0.1/build -DMPI_DIR=/usr/local/mpi3 
-DP4EST_DIR=../p4est/p4est-2.0 -DPETSC_DIR=../new-petsc-3.7.7/petsc-3.7.7 
-DHDF5_DIR=/usr/lib64/ -DMETIS_DIR=../metis-5.1.0/metis-5.1.0 
-DDEAL_II_WITH_METIS=ON -DDEAL_II_WITH_HDF5=OFF -DDEAL_II_WITH_PETSC=ON 
-DDEAL_II_WITH_P4EST=OFF -DDEAL_II_WITH_MPI=ON 
-DDEAL_II_FORCE_BUNDLED_BOOST=ON-DDEAL_II_ALLOW_AUTODETECTION=ON

$make

$make install

But after "make install", the installation was interupped, and the errors 
are showed as:

Scanning dependencies of target obj_rol_debug
[ 82%] Built target obj_rol_debug
Scanning dependencies of target step-3.release
[ 82%] Building CXX object 
examples/CMakeFiles/step-3.release.dir/step-3/step-3.cc.o
[ 82%] Linking CXX executable ../bin/step-3.release
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 1>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<3, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<3, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 1>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<1, 2>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 3>::StraightBoundary()'
../lib/libdeal_II.so.9.0.1: undefined reference to 
`dealii::StraightBoundary<2, 2>::StraightBoundary()'
make[2]: *** [bin/step-3.release] Error 1
make[1]: *** [examples/CMakeFiles/step-3.release.dir/all] Error 2
make: *** [all] Error 2

So, I don't what happened and know how to do. Hope some suggestions! Thanks 
in advance!

Best,
Chucui




-- 
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: How to plot a function including "arctan" correctly?

[chucui1016]

Dear Wolfgang,

Thank you very much!


>> So, my questions are: 
>> 1.If I need to get 
>>
>> 4.JPG 
>> 
>> 
>>in every time step, which means that phi_x and phi_y are updated in 
>> every time step, and the information of phi_x and phi_y are calculated 
>> from "phi" by using function: "get_function_gradients". 
>> 
>> And I don't know the exact value of phi and grad_phi in every time step. 
>> So how can I remove  oscillations? 

> I don't understand this. If you don't know what grad phi is, how do you 
> want to compute atan2(phi_y,phi_x) and from that F? If you don't know 
> what phi is, you also can't compute grad phi via a finite different 
> approximation. 

> It doesn't matter whether you compute the gradient exactly or via a FDM 
> approximation *conceptually*, but it will almost certainly be *more 
> accurate* to use the exact derivative with get_function_gradients(). 

I want to describe what do I need the gradient for: 
https://groups.google.com/forum/#!topic/dealii/YSzg4k47T00
As you can see, I want to update \phi in every time step, then use 
"get_function_gradients()" to get grad_phi in every time step. But my 
problem is: 
In 0th time step, I have exact original phi function, then project it onto 
the C0 finite element space and get the degree of freedom vector: phi_0. 
Then I use "get_function_gradients()" to get grad_phi, but the result is 
not so good as above, unless use FDM to approximate grad_phi.
For these experiences in 0th time step, I guess that maybe I need to do 
something like this in every time step. I get that 
"get_function_gradients()" is more accurate, but it can't deal with some 
oscillations I think. As a result, I want to use FDM  to approximate 
grad_phi in every time step. So I need the exact real coordinates of grid 
nodes (like exact (x,y), not only the number of the cell and the number 
which point of quadrature_point), then get the real coordinates of adjacent 
nodes, and get the values of phi_up, phi_down, phi_left, phi_right.
Now I use 
  cell_matrix(i,j) += (fe_values.shape_value (i, q_index) *
   fe_values.shape_value (j, q_index) *
   fe_values.JxW (q_index));
to describe every point I can get their value, how can I get the exact real 
coordinates of thses points (Gauss points) or grid nodes? Which function 
can return the exact real coordinates?

Thank you very much!

Best, 
Chucui

-- 
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: Behavior of FESystem with 2 vectors and 2 scalars

[chucui1016]

Dear Daniel，

Thank you very much! My questions are all solved. 

Best,
Chucui

-- 
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: Behavior of FESystem with 2 vectors and 2 scalars

[chucui1016]

The attachment are 2 code files：LittleCase0-4.cc (for test), 
LittleCase0-5.cc (RT element).

在 2018年8月26日星期日 UTC+8下午8:32:40，chucu...@gmail.com写道：
>
> Hi, All,
>
> I have a question of behavior of FESystem with 2 vectors and 2 scalars, my 
> original problem is to solve:
>
> [image: 56.JPG]
> for \phi, \sigma, \mu, \gamma, I choose FESystem to describe the 
> information of them
> FESystemfe;
> [...]
> fe (FE_DGQ(degree), 1,
> FE_RaviartThomas(degree), 1,
> FE_DGQ(degree), 1,
> FE_RaviartThomas(degree), 1),
>
> The test code is LittleCase0-5.cc, and after make run, the error 
> information is: 
> Linking CXX executable Case0-5
> [ 50%] Built target Case0-5
> [100%] Run with Debug configuration
> before_start0
> dof 9-1 
> dof 9-2 6|112
> dof 9-3 
> dof 9-4 
> dof 9-5 
>
> 
> An error occurred in line <1746> of file 
>  in function
> void dealii::DoFTools::count_dofs_per_block(const DoFHandlerType&, 
> std::vector&, const std::vector&) [with 
> DoFHandlerType = dealii::DoFHandler<2>]
> The violated condition was: 
> target_block.size()==fe.n_blocks()
> Additional information: 
> Dimension 6 not equal to 4.
>
>
> to test the behavior of FESystem with 2 vectors and 2 scalars, I choose 
> another FESystem:
> FESystemfe;
> [...]
> fe (FE_Q(degree), 1,
> FE_Q(degree), dim,
> FE_Q(degree), 1,
> FE_Q(degree), dim),
>
> the code is LittleCase0-4.cc, and the result after "make run" is:
> Linking CXX executable Case0-4
> [ 50%] Built target Case0-4
> [100%] Run with Debug configuration
> before_start0
> dof 9-1 
> dof 9-2 6|54
> dof 9-3 
> dof 9-4 
> dof 9-5 
> dof 9-6 
> dof 9-7 
>Number of active cells: 4
>Number of degrees of freedom: 54 (9+27|9+9)
> dof 9-8 
> dof 9-9 
> dof 9-10 
> dof 9-11 
> dof 9-12 
> dof 9-13 
> setup_dofs0
> [100%] Built target run
>
>
> My question is:
>
> 1. Can FESystem countain 2 vectors and 2 scalars? If so, how to make 
> "dof_handler", "dofs_per_block", and "block_component (target_block)" 
> correctly matching? As above, I make some faults. For example, in 
> LittleCase0-4.cc, I use this defination for "block_component" 
> (target_block):
> std::vector block_component (2*dim+2,1);
> block_component[0] = 0;
> block_component[1] = 1;
> block_component[dim+1] = 2;//=3-1
> block_component[dim+2] = 3;
> and te result is
> Number of degrees of freedom: 54 (9+27|9+9)
> but it should be: " 54 (9+18 | 9+18) ", I don't know how to control it as 
> I cannot set  "std::vector block_component (2*dim+2,1,3);", 
> where "1" and "3" means that the vector block should be in the first and 
> third block in the whole.
>
> 2. Can I use RT element by this way (in Block)? If so, how to correct the 
> LittleCase0-5.cc? If not, how to implement the original problem? Maybe use 
> 2 FESystem fe_1, fe_2, both of which contain only 1 vector (RT element) and 
> 1 scalar (DG element)?
>
>
> Thanks in advance!
>
> Best,
> Chucui
>
>
>
>
>  
>
>
>
>
>

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

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

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

#include 
#include 
#include 
#include 

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

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include 

// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include 
#include
#include "stdio.h"
#include
// This is C++:
#include 
#include 
#include 
#define random(x) (rand()%x)
#define MAX 2
// As in all programs, the namespace dealii is included:
namespace Step22
{
  using namespace dealii;


  template 
  class StokesProblem
  {
  public:
StokesProblem (const unsigned int degree);
void run ();

  private:
void setup_dofs ();

void refine_mesh ();



const unsigned int   degree;

Triangulation   triangulation_active;
Triangulation   triangulation;
FESystemfe;
  

[deal.II] Behavior of FESystem with 2 vectors and 2 scalars

[chucui1016]

Hi, All,

I have a question of behavior of FESystem with 2 vectors and 2 scalars, my 
original problem is to solve:

[image: 56.JPG] 
for \phi, \sigma, \mu, \gamma, I choose FESystem to describe the 
information of them
FESystemfe;
[...]
fe (FE_DGQ(degree), 1,
FE_RaviartThomas(degree), 1,
FE_DGQ(degree), 1,
FE_RaviartThomas(degree), 1),

The test code is LittleCase0-5.cc, and after make run, the error 
information is: 
Linking CXX executable Case0-5
[ 50%] Built target Case0-5
[100%] Run with Debug configuration
before_start0
dof 9-1 
dof 9-2 6|112
dof 9-3 
dof 9-4 
dof 9-5 


An error occurred in line <1746> of file 
 in function
void dealii::DoFTools::count_dofs_per_block(const DoFHandlerType&, 
std::vector&, const std::vector&) [with 
DoFHandlerType = dealii::DoFHandler<2>]
The violated condition was: 
target_block.size()==fe.n_blocks()
Additional information: 
Dimension 6 not equal to 4.


to test the behavior of FESystem with 2 vectors and 2 scalars, I choose 
another FESystem:
FESystemfe;
[...]
fe (FE_Q(degree), 1,
FE_Q(degree), dim,
FE_Q(degree), 1,
FE_Q(degree), dim),

the code is LittleCase0-4.cc, and the result after "make run" is:
Linking CXX executable Case0-4
[ 50%] Built target Case0-4
[100%] Run with Debug configuration
before_start0
dof 9-1 
dof 9-2 6|54
dof 9-3 
dof 9-4 
dof 9-5 
dof 9-6 
dof 9-7 
   Number of active cells: 4
   Number of degrees of freedom: 54 (9+27|9+9)
dof 9-8 
dof 9-9 
dof 9-10 
dof 9-11 
dof 9-12 
dof 9-13 
setup_dofs0
[100%] Built target run


My question is:

1. Can FESystem countain 2 vectors and 2 scalars? If so, how to make 
"dof_handler", "dofs_per_block", and "block_component (target_block)" 
correctly matching? As above, I make some faults. For example, in 
LittleCase0-4.cc, I use this defination for "block_component" 
(target_block):
std::vector block_component (2*dim+2,1);
block_component[0] = 0;
block_component[1] = 1;
block_component[dim+1] = 2;//=3-1
block_component[dim+2] = 3;
and te result is
Number of degrees of freedom: 54 (9+27|9+9)
but it should be: " 54 (9+18 | 9+18) ", I don't know how to control it as I 
cannot set  "std::vector block_component (2*dim+2,1,3);", 
where "1" and "3" means that the vector block should be in the first and 
third block in the whole.

2. Can I use RT element by this way (in Block)? If so, how to correct the 
LittleCase0-5.cc? If not, how to implement the original problem? Maybe use 
2 FESystem fe_1, fe_2, both of which contain only 1 vector (RT element) and 
1 scalar (DG element)?


Thanks in advance!

Best,
Chucui




 




-- 
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] If I define two FESystem , how to use "FEValuesExtractors" separately？

[chucui1016]

Hi, all

I have a simple question: if I define two FESystem:
[...]
FESystemfe_1;
FESystemfe_2;
[...]
fe_1 (FE_Q(degree), 1,
  FE_Q(degree), 1),
fe_2 (FE_RaviartThomas(degree), 1),
[...]

when I need to use  "FEValuesExtractors" separately, how to define which 
FEValues are subscripts in operator[] expressions on? fe_1 or fe_2?
if my code 
FEValuesExtractors::Scalar fe_1_0(0);
FEValuesExtractors::Scalar fe_1_1(1);

FEValuesExtractors::Vector fe_2_vector(0);

can "FEValuesExtractors" match to the corresponding FESystem （fe_1 or 
fe_2）automatically?

Thanks very much!

Best,
Chucui

-- 
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: implemention of periodic boundary condition

[chucui1016]

Dear Daniel and Wolfgang,

Thank you very much! In fact, I only want the gradient of periodic, not the 
Dirichlet periodic boundary. I thought they are same in periodic boundary 
conditions before, i.e. all informations of the solution on the face pairs 
are same. Now, I get that I have made mistake.

Why I need the gradient of periodic only? Because in my case, I have some 
boundary integrals in my cost functional, I need to make all of them equal 
to zero, so I can choose the gradient of periodic only.

But in the analysis of weak solution of original PDEs in my case, the space 
of trial and test funcution are (H^1) * (H^1) * (L^2) * (H^1) (sobolev 
space), as there are the gradient of 1st, 2nd, 4th component in the weak 
form. And I choose C0 * C0 * C0 * C0 FEM space

So I cannot make sure the continuity of the gradient in boundary or 
internal faces. But I don' know how to implement what I need? Maybe some 
penalty part like in Discontinuous Gradient Method to make up for 
continuity of the gradient?

Thank you very much!

Best,
Chucui

-- 
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: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

Thank you very much! I can test my period boundary condition added truly. 
Your advision helps me a lot !

As my question is about periodic boundary condition like 

[image: 42.JPG] 
just like Dirichlet boundary condition.

But if "periodic boundary condition" means that all information of solution 
on periodic faces pair are equal, not only the value of solution (like 
Dirichlet boundary condition above), but also gradient value of solution 
(like Neuman boundary condition, as normal vector is opposite, so if 
grad_u_on_left = grad_u_on_right, for example, so n_left * grad_u_on_left + 
n_right * grad_u_on_right, top and bottom as well) . In my case, periodic 
boundary condition should include

[image: 46.JPG] 



Does deal.ii implement this kind of periodic boundary condition in fact? Is 
it different from what we have done above for periodic boundary condition 
like Dirichlet?

Thank you very much!

Best,
Chucui


-- 
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: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

Thank you very much! I use code like this:
constraints.distribute_local_to_global(local_matrix, local_dof_indices, 
system_matrix) as you say, 

and use codes like this
constraints.distribute(solution);
after every "solve()" step.

and now, my LittleCase0-3.cc can run. 
Maybe a question more : how to make sure your period boundary coditions are 
added truly?

Thank you very much!

Best, 
Chucui

-- 
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: implemention of periodic boundary condition

[chucui1016]

Dear Jean-Paul,

Thank you very much for your patience, when I print some information like 
this:

template 
void Problem::run ()
{
[...]
  setup_dofs ();
std::cout << "test-1 "<< std::endl;  
VectorTools::project (dof_handler_0, constraints_phi, 
QGauss(degree+2),
SolutionsPhi(),
phi_0);  
std::cout << "test-2 " << std::endl;   
   
VectorTools::project (dof_handler_2, constraints_q, 
QGauss(degree+2),
SolutionsQ(),
q_0);  
 std::cout << "test-3 " << std::endl;   
assemble_system_1 (phi_0);
 std::cout << "test-4" << std::endl;   
assemble_rhs_1 (phi_0, q_0);
 std::cout << "test-5" << std::endl; 

and the output is

Linking CXX executable Case0-3
[ 50%] Built target Case0-3
[100%] Run with Debug configuration
   Number of active cells: 16384
   Number of degrees of freedom: 98818 (16641+16641+65536)
test-1 
test-2 
test-3 


An error occurred in line <1764> of file 
 
in function
void 
dealii::SparseMatrix::add(dealii::SparseMatrix::size_type, 
dealii::SparseMatrix::size_type, number) [with number = double; 
dealii::SparseMatrix::size_type = unsigned int]
The violated condition was: 
(index != SparsityPattern::invalid_entry) || (value == number())
Additional information: 
You are trying to access the matrix entry with index <0,1>, but this 
entry does not exist in the sparsity pattern of this matrix.


So it can't run when arrive at "assemble_system_1", so "setup-dof()", 
 "assemble_system_1()", "run()" are left—— this is LittleCase0-3.cc
if only "setup_dof()", "run()" are left——this is LittleCase0-4.cc

Hope this is a minimal example. Thank you very much!

Best,
Chucui


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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include 

// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include 

// This is C++:
#include 
#include 
#include 

// As in all programs, the namespace dealii is included:
namespace Step22
{
  using namespace dealii;


  template 
  class Problem
  {
  public:
Problem (const unsigned int degree);
void run ();

  private:
void setup_dofs ();

void assemble_system_1 (const Vector phi_0);

const unsigned int   degree;

Triangulation   triangulation;
Triangulation   triangulation_active;

FESystemfe;
FE_Qfe_phi;
FE_Qfe_mu;
FE_DGQ  fe_q;
DoFHandler  dof_handler;
DoFHandler  dof_handler_0;
DoFHandler  dof_handler_1;
DoFHandler  dof_handler_2;

ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;

BlockSparsityPattern  sparsity_pattern;


BlockSparseMatrix system_matrix;
BlockSparseMatrix system_matrix_1;
BlockSparseMatrix free_energy_matrix;

Vector phi_0, q_0;

BlockVector solution;
BlockVector old_solution;
BlockVector older_solution;
BlockVector system_rhs;
BlockVector system_rhs_1;
BlockVector solution_1;

double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;

ConvergenceTable convergence_table;

  };
  

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

virtual double value (const Point   &

[deal.II] Re: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

I'm sorry to make this mistake, (and thank Jean-Paul very much). The newest 
code focus on build the sparsity pattern in "setup_dof()", while other 
parts (assemble_system(), assemble_rhs(), run() ) need to be left to make 
the code work and maintain the structure of 3 components of the blockvector 
solution.

Thank you much!

Best,
Chucui

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include 

// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include 

// This is C++:
#include 
#include 
#include 

// As in all programs, the namespace dealii is included:
namespace Step22
{
  using namespace dealii;


  template 
  class Problem
  {
  public:
Problem (const unsigned int degree);
void run ();

  private:
void setup_dofs ();
void assemble_system (const BlockVector old_solution,
  const BlockVector older_solution);
void assemble_rhs (const double time,
   const BlockVector old_solution,
   const BlockVector older_solution);
void solve ();

double free_energy_compute (const BlockVector solution);

void assemble_system_1 (const Vector phi_0);
void assemble_rhs_1 (const Vector phi_0, const Vector q_0);
void solve_1 ();

const unsigned int   degree;

Triangulation   triangulation;
Triangulation   triangulation_active;

FESystemfe;
FE_Qfe_phi;
FE_Qfe_mu;
FE_DGQ  fe_q;
DoFHandler  dof_handler;
DoFHandler  dof_handler_0;
DoFHandler  dof_handler_1;
DoFHandler  dof_handler_2;

ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;

BlockSparsityPattern  sparsity_pattern;


BlockSparseMatrix system_matrix;
BlockSparseMatrix system_matrix_1;
BlockSparseMatrix free_energy_matrix;

Vector phi_0, q_0;

BlockVector solution;
BlockVector old_solution;
BlockVector older_solution;
BlockVector system_rhs;
BlockVector system_rhs_1;
BlockVector solution_1;

double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;

ConvergenceTable convergence_table;

  };
  

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

virtual double value (const Point   &p,
  const unsigned int  component = 0) const;
virtual Tensor<1,dim> gradient (const Point   &p,
const unsigned int  component = 0) const;
  };
  
  
  template 
  double SolutionsPhi::value (const Point   &p,
   const unsigned int component) const
  {
 //if (component == 0)
 double phi_value = 0;

  Tensor<1,dim> r;

  r[0] = std::sqrt((p[0] - 0.3 + 0.01) * (p[0] - 0.3 + 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
  r[1] = std::sqrt((p[0] - 0.7 - 0.01) * (p[0] - 0.7 - 0.01) + (p[1] - 0.5) * (p[1] - 0.5)) ;
  
  if ((r[0] >= (0.2 + 0.01)) && (r[1] >= (0.2 + 0.01)))
   {phi_value = -1;
   }
  else if (r[0] < (0.2 + 0.01))
{phi_value = std::tanh((0.2-std::sqrt((p[0]-0.3+0.01) * (p[0]-0.3+0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
}
   else if (r[1] < (0.2 + 0.01))
 {phi_value = std::tanh((0.2-std::sqrt((p[0]-0.7-0.01) * (p[0]-0.7-0.01) + (p[1] -0.5) * (p[1] -0.5)))/0.01);
 } 
else
   

[deal.II] Re: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

Thank you very much! The LittleCase0-3.cc is my work case with the same 
problem as I say.

Simply describe what this code work for: to solve



from



where the three variables have periodic boundary conditions on the domain 
[0,1]*[0,1]. So the weak form has no part of boundary integral.

Thank you very much!

Best,
Chucui

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include 

// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include 

// This is C++:
#include 
#include 
#include 

// As in all programs, the namespace dealii is included:
namespace Step22
{
  using namespace dealii;


  template 
  class Problem
  {
  public:
Problem (const unsigned int degree);
void run ();

  private:
void setup_dofs ();
void assemble_system (const BlockVector old_solution,
  const BlockVector older_solution);
void assemble_rhs (const double time,
   const BlockVector old_solution,
   const BlockVector older_solution);
void solve ();
void output_results_in_routine (const unsigned int timestep_number)  const;
void output_results (const unsigned int refinement_cycle) const;
void output_results_phi ()  const;
void output_results_q ()  const;
void output_results_initial ()  const;

void refine_mesh ();
void process_solution(const double time);
double free_energy_compute (const BlockVector solution);

void assemble_system_1 (const Vector phi_0);
void assemble_rhs_1 (const Vector phi_0, const Vector q_0);
void solve_1 ();

const unsigned int   degree;

Triangulation   triangulation;
Triangulation   triangulation_active;

FESystemfe;
FE_Qfe_phi;
FE_Qfe_mu;
FE_DGQ  fe_q;
DoFHandler  dof_handler;
DoFHandler  dof_handler_0;
DoFHandler  dof_handler_1;
DoFHandler  dof_handler_2;

ConstraintMatrix constraints;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_mu;
ConstraintMatrix constraints_q;

BlockSparsityPattern  sparsity_pattern;
//BlockSparsityPattern  sparsity_pattern_mu;

BlockSparseMatrix system_matrix;
BlockSparseMatrix system_matrix_1;
BlockSparseMatrix free_energy_matrix;

//SparseMatrix system_matrix_mu_0;
Vector phi_0, q_0;

BlockVector solution;
BlockVector old_solution;
BlockVector older_solution;
BlockVector system_rhs;
BlockVector system_rhs_1;
BlockVector solution_1;

double time_step, time;
const double D, eps, lambda;//, s;
unsigned int timestep_number;

ConvergenceTable convergence_table;

  };
  

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

virtual double value (const Point   &p,
  const unsigned int  component = 0) const;
virtual Tensor<1,dim> gradient (const Point   &p,
const unsigned int  component = 0) const;
  };
  
  
  template 
  double SolutionsPhi::value (const Point   &p,
   const unsigned int component) const
  {
 //if (component == 0)
 double phi_value = 0;
// phi_value = std::tanh((0.2-std::

[deal.II] Re: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

Thank you very much for your rapid reply! I use overload 

 as 
you say, and it can run before assemble_system, the relevant part of my 
code is
  template 
  void Problem::run ()
  {
   
const unsigned int n_cycles = 1;
double entropy = 0;

for (unsigned int cycle=0; cycle::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f::faces_per_cell; 
++f)
{
   if (cell->face(f)->at_boundary())
   {  
  if (std::fabs(cell->face(f)->center()(0) - 
(2.25)) < 1e-12)
  
  cell->face(f)->set_boundary_id (1);
  else if(std::fabs(cell->face(f)->center()(1) - 
(0)) < 1e-12)
   
   cell->face(f)->set_boundary_id (2);
   else if(std::fabs(cell->face(f)->center()(1) 
- (2.25)) < 1e-12)
 
 cell->face(f)->set_boundary_id (3);
 else 
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
  {
  cell->face(f)->set_boundary_id 
(0);
  }  
   } 
   }
}
 
 }
else
{
//triangulation.refine_global (1);
}
   setup_dofs ();
   std::cout << "start" << cycle << std::endl;  

   
VectorTools::project (dof_handler_phi, constraints, 
QGauss(degree+2),
InitialValues_phi(),
phi_0);   
VectorTools::project (dof_handler_e, constraints, 
QGauss(degree+2),
InitialValues_e(),
e_0);  
std::cout << "phi_0" << cycle << std::endl;
//VectorTools::project (dof_handler_U, constraints_U, 
QGauss(degree+2),
//InitialValues_U(),
//U_n);
//std::cout << "U_0" << cycle << std::endl;  
VectorTools::project (dof_handler_q, constraints_q, 
QGauss(degree+2),
InitialValues_q(),
q_0);  
std::cout << "q_0" << cycle << std::endl; 
assemble_system_1 (phi_0, e_0);
[...]
and my setup_dofs ():
  template 
  void Problem::setup_dofs ()
  {
std::vector block_component (4);
block_component[0] = 0;
block_component[1] = 1; 
block_component[2] = 2;
block_component[3] = 3;

dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, 
component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, 
component_mask);

constraints.close (); 

std::vector dofs_per_block (4); 
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, 
block_component);

const unsigned int n_phi = dofs_per_block[0],
   n_e = dofs_per_block[1],
   n_U = dofs_per_block[2],
   n_q = dofs_per_block[3];
system_matrix_1.clear ();
   
BlockDynamicSparsityPattern dsp (4,4);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);  
dsp.block(2,0).reinit (n_U, n_phi);
dsp.block(3,0).reinit (n_q, n_phi); 
 
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
dsp.block(2,1).reinit (n_U, n_e);
dsp.block(3,1).reinit (n_q, n_e); 

dsp.block(0,2).reinit (n_phi, n_U);
dsp.block(1,2).reinit (n_e, n_U);  
dsp.block(2,2).reinit (n_U, n_U);
dsp.block(3,2).reinit (n_q, n_U); 
 
dsp.block(0,3).reinit (n_phi, n_q);
dsp.block(1,3).reinit (n_e, n_q);
dsp.block(2,3).reinit (n_U, n_q);
dsp.block(3,3).reinit (n_q, n_q); 

dsp.collect_sizes();

DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);

system_matrix_1.reinit (sparsity_pattern);
[...]
the output and error is
start0
phi_00
q_00


An error occurred in line <1764> of file 
 in 
function
void 
dealii::SparseMatrix::add(dealii::SparseMatrix::size_type, 
dealii::SparseMatrix::size_type, number) [with

[deal.II] Re: implemention of periodic boundary condition

[chucui1016]

Dear Daniel,

Thank you very much for your rapid reply! I use overload 

 as 
you say, and it can run before assemble_system, the relevant part of my 
code is
  template 
  void Problem::run ()
  {
   
const unsigned int n_cycles = 1;
double entropy = 0;

for (unsigned int cycle=0; cycle::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f::faces_per_cell; 
++f)
{
   if (cell->face(f)->at_boundary())
   {  
  if (std::fabs(cell->face(f)->center()(0) - 
(2.25)) < 1e-12)
  
  cell->face(f)->set_boundary_id (1);
  else if(std::fabs(cell->face(f)->center()(1) - 
(0)) < 1e-12)
   
   cell->face(f)->set_boundary_id (2);
   else if(std::fabs(cell->face(f)->center()(1) 
- (2.25)) < 1e-12)
 
 cell->face(f)->set_boundary_id (3);
 else 
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
  {
  cell->face(f)->set_boundary_id 
(0);
  }  
   } 
   }
}
 
 }
else
{
//triangulation.refine_global (1);
}
   setup_dofs ();
   std::cout << "start" << cycle << std::endl;  

   
VectorTools::project (dof_handler_phi, constraints, 
QGauss(degree+2),
InitialValues_phi(),
phi_0);   
VectorTools::project (dof_handler_e, constraints, 
QGauss(degree+2),
InitialValues_e(),
e_0);  
std::cout << "phi_0" << cycle << std::endl;
//VectorTools::project (dof_handler_U, constraints_U, 
QGauss(degree+2),
//InitialValues_U(),
//U_n);
//std::cout << "U_0" << cycle << std::endl;  
VectorTools::project (dof_handler_q, constraints_q, 
QGauss(degree+2),
InitialValues_q(),
q_0);  
std::cout << "q_0" << cycle << std::endl; 
assemble_system_1 (phi_0, e_0);
[...]
and my setup_dofs ():
  template 
  void StokesProblem::setup_dofs ()
  {
std::vector block_component (4);
block_component[0] = 0;
block_component[1] = 1; 
block_component[2] = 2;
block_component[3] = 3;

dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, 
component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, 
component_mask);

constraints.close (); 

std::vector dofs_per_block (4); 
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, 
block_component);

const unsigned int n_phi = dofs_per_block[0],
   n_e = dofs_per_block[1],
   n_U = dofs_per_block[2],
   n_q = dofs_per_block[3];
system_matrix_1.clear ();
   
BlockDynamicSparsityPattern dsp (4,4);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);  
dsp.block(2,0).reinit (n_U, n_phi);
dsp.block(3,0).reinit (n_q, n_phi); 
 
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
dsp.block(2,1).reinit (n_U, n_e);
dsp.block(3,1).reinit (n_q, n_e); 

dsp.block(0,2).reinit (n_phi, n_U);
dsp.block(1,2).reinit (n_e, n_U);  
dsp.block(2,2).reinit (n_U, n_U);
dsp.block(3,2).reinit (n_q, n_U); 
 
dsp.block(0,3).reinit (n_phi, n_q);
dsp.block(1,3).reinit (n_e, n_q);
dsp.block(2,3).reinit (n_U, n_q);
dsp.block(3,3).reinit (n_q, n_q); 

dsp.collect_sizes();

DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);

system_matrix_1.reinit (sparsity_pattern);
[...]
the output and error is
start0
phi_00
q_00


An error occurred in line <1764> of file 
 in 
function
void 
dealii::SparseMatrix::add(dealii::SparseMatrix::size_type, 
dealii::SparseMatrix::size_type, number)

[deal.II] Re: implemention of periodic boundary condition

[chucui1016]

Dear Jean-Paul,

Thank you very much! I had read step-45, but still have some questions:
1. How to use function "DoFTools::make_periodicity_constraints() "  only? 
if my code need not to compute parallel case.  If I choose
void DoFTools::make_periodicity_constraints ( const FaceIterator & face_1,
const typename identity 
< 
FaceIterator >::type & face_2,
AffineConstraints 
< 
double > & constraints,
const ComponentMask 
 & 
component_mask 
= ComponentMask 
(),
const bool face_orientation = true,
const bool face_flip = false,
const bool face_rotation = false,
const FullMatrix 
< double 
> & matrix = FullMatrix 

(),
const std::vector< unsigned int > & first_vector_components = std::vector<
unsigned int>() 
)
how to choose face1 and face2 we need there correctly? If I use 
template
void DoFTools::make_periodicity_constraints ( const std::vector< GridTools::
PeriodicFacePair 

< typename DoFHandlerType::cell_iterator >> & periodic_faces,
AffineConstraints 
< 
double > & constraints,
const ComponentMask 
 & 
component_mask 
= ComponentMask 
(),
const std::vector< unsigned int > & first_vector_components = std::vector<
unsigned int>() 
)
how to make a 
const std::vector< GridTools::PeriodicFacePair 

< typename DoFHandlerType::cell_iterator >> & periodic_faces
maybe I still need to use "GridTools::collect_periodic_faces" to input 
information of faces into PeriodicFacePair. If I do this,  I am not use 
"DoFTools::make_periodicity_constraints() "  only


2. In step-45, the complete implemention is parallel case. For every step, 
it use functions used in parallel case. If I want to implement a case not 
in parallel, shall I rewrite all of them and remove the parallel parts 
(still not succeed), or use it directly, but in shell run my program via
mpirun -np 1 ./myCase


Thank you very much!

Best,
Chucui

-- 
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: implemention of periodic boundary condition

[chucui1016]

Hi All,

I am working on this project and find some problems new:

1. Why the source code of function "GridTools::collect_periodic_faces" use "
cell_iterator" rather than "active_cell_iterator"? As I need to refine my 
domain by 
triangulation.refine_global (7);
in https://www.mail-archive.com/dealii@googlegroups.com/msg00220.html
I find the same problem, what "Periodic boundaries can only have a 
difference of 1 refinement level between pairs of faces." means?  And how 
can I implement my project for refine global (not in parallel)?

2. If I use "Triangulation" rather than 
"parallel::distributed::Triangulation", shall I use the function 
"GridTools::collect_periodic_faces"? Or just use only 
"DoFTools::make_periodicity_constraints"?

3. If use "DoFTools::make_periodicity_constraints" only, my code is
  template 
  void StokesProblem::setup_dofs ()
  {

std::vector block_component (4);
block_component[0] = 0;
block_component[1] = 1; 
block_component[2] = 2;
block_component[3] = 3;


dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);

 
constraints.clear ();
const unsigned int n_dofs_per_face = fe.dofs_per_face;
FullMatrix rotation_matrix(n_dofs_per_face);
for (unsigned int d=0; d offset;
const ComponentMask component_mask = ComponentMask();
   
   std::vector first_vector_components;
   first_vector_components.push_back(0);
   first_vector_components.push_back(1);
   first_vector_components.push_back(2);
   first_vector_components.push_back(3);
   
 const typename DoFHandler::active_face_iterator face_0; 
 const typename DoFHandler::active_face_iterator face_1;
 const typename DoFHandler::active_face_iterator face_2;
 const typename DoFHandler::active_face_iterator face_3;

 for (typename DoFHandler::active_cell_iterator cell = 
dof_handler.begin_active();
  cell != dof_handler.end();
  ++cell)
   { 

 for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; 
++i)
   {
 const typename DoFHandler::active_face_iterator face = 
cell->face(i);
 
 if (face->at_boundary() && face->boundary_id() == 0)
   {
const typename DoFHandler::active_face_iterator face_0 
= face; 
   }
 
 else if (face->at_boundary() && face->boundary_id() == 1)
   {
 const typename DoFHandler::active_face_iterator 
face_1 = face; 
   }
 else if (face->at_boundary() && face->boundary_id() == 2)
   {
const typename DoFHandler::active_face_iterator face_2 
= face; 
   }
 
else if (face->at_boundary() && face->boundary_id() == 3)
   {
 const typename DoFHandler::active_face_iterator 
face_3 = face; 
   }   
   
   }
 
 DoFTools::make_periodicity_constraints (face_0, face_1, 
constraints, component_mask, true, false, false, rotation_matrix, 
first_vector_components);
 DoFTools::make_periodicity_constraints (face_2, face_3, 
constraints, component_mask, true, false, false, rotation_matrix, 
first_vector_components);
   }
[...]

but get the error:
Linking CXX executable Problem13-2
/.../13-4-pbc-for-run/together-pbc.cc:647: error: undefined reference to 
'void 
dealii::DoFTools::make_periodicity_constraints, false> > 
>(dealii::TriaActiveIterator, false> > const&, 
dealii::identity, false> > >::type const&, 
dealii::ConstraintMatrix&, dealii::ComponentMask const&, bool, bool, bool, 
dealii::FullMatrix const&, std::vector > const&)'
/.../13-4-pbc-for-run/together-pbc.cc:648: error: undefined reference to 
'void 
dealii::DoFTools::make_periodicity_constraints, false> > 
>(dealii::TriaActiveIterator, false> > const&, 
dealii::identity, false> > >::type const&, 
dealii::ConstraintMatrix&, dealii::ComponentMask const&, bool, bool, bool, 
dealii::FullMatrix const&, std::vector > const&)'
collect2: error: ld returned 1 exit status
make[3]: *** [Problem13-2] Error 1
make[2]: *** [CMakeFiles/Problem13-2.dir/all] Error 2
make[1]: *** [CMakeFiles/run.dir/rule] Error 2
make: *** [run] Error 2


I wander what is the reason for this issue.

Thanks very much!

Best,
Chucui

在 2018年8月5日星期日 UTC+8下午3:28:01，chucu...@gmail.com写道：
>
> Hi All,
>
> I am having a problem of implementing periodic boundary condition, I have 
> 4 variables to solve, and I make them into a blockvector:
>
>
> 
>
> each component of it is a scalar. All of them have periodic boundary 
> condition: 
>
>
> 

[deal.II] implemention of periodic boundary condition

[chucui1016]

Hi All,

I am having a problem of implementing periodic boundary condition, I have 4 
variables to solve, and I make them into a blockvector:



each component of it is a scalar. All of them have periodic boundary 
condition: 



on domain:



and my code is:
  template 
  void Problem::setup_dofs ()
  {

std::vector block_component (4);
block_component[0] = 0;
block_component[1] = 1; 
block_component[2] = 2;
block_component[3] = 3;


dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);

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

  Tensor<1,dim> offset;

  std::vector::cell_iterator> >
  periodicity_vector;

const unsigned int direction = 0;

   const ComponentMask component_mask = ComponentMask();
   
   std::vector first_vector_components;
   first_vector_components.push_back(0);
   first_vector_components.push_back(1);
   first_vector_components.push_back(2);
   first_vector_components.push_back(3);
   
   GridTools::collect_periodic_faces(dof_handler,
  /*b_id1*/ 0,
  /*b_id2*/ 1,
  /*direction*/ 0,
  periodicity_vector);
   GridTools::collect_periodic_faces(dof_handler,
  /*b_id1*/ 2,
  /*b_id2*/ 3,
  /*direction*/ 0,
  periodicity_vector);
   DoFTools::make_periodicity_constraints >
  (periodicity_vector, constraints);

constraints.close ();

[...]  

  template 
  void Problem::run ()
  {
   
const unsigned int n_cycles = 1;

  for (unsigned int cycle=0; cycle::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
for (unsigned int f=0; f::faces_per_cell; 
++f)
   if (cell->face(f)->at_boundary())

  if (std::fabs(cell->face(f)->center()(0) - 
(2.25)) < 1e-12)
  cell->face(f)->set_boundary_id (1);
  else if(std::fabs(cell->face(f)->center()(1) - 
(0)) < 1e-12)
   cell->face(f)->set_boundary_id (2);
else 
if(std::fabs(cell->face(f)->center()(1) - (2.25)) < 1e-12)
 cell->face(f)->set_boundary_id (3);
 else 
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
  cell->face(f)->set_boundary_id 
(0);

   std::vector<
   GridTools::PeriodicFacePair::cell_iterator>>
   periodicity_vector;
FullMatrix rotation_matrix(dim);
rotation_matrix[0][0]=1.;
rotation_matrix[1][1]=1.;
   
   const ComponentMask component_mask = ComponentMask();
   
   std::vector first_vector_components;
   first_vector_components.push_back(0);
   first_vector_components.push_back(1);
   first_vector_components.push_back(2);
   first_vector_components.push_back(3);
  GridTools::collect_periodic_faces(dof_handler,
  /*b_id1*/ 0,
  /*b_id2*/ 1,
  /*direction*/ 0,
  periodicity_vector);
  GridTools::collect_periodic_faces(dof_handler,
  /*b_id1*/ 2,
  /*b_id2*/ 3,
  /*direction*/ 0,
  periodicity_vector); 
   DoFTools::make_periodicity_constraints >
  (periodicity_vector, constraints);
 }
}

setup_dofs ();
[...]
but when make run it, got the following error:
An error occurred in line <3751> of file 
<.../boost/deal.ii-8.5.1/src/source/grid/grid_tools.cc> in function
void dealii::GridTools::collect_periodic_faces(const MeshType&, 
dealii::types::boundary_id, dealii::types::boundary_id, int, 
std::vector >&, const dealii::Tensor<1, MeshType:: 
space_dimension>&, const dealii::FullMatrix&) [with MeshType = 
dealii::DoFHandler<2>; dealii::types::boundary

[deal.II] Re: when make, something about libgfortran may be wrong

[chucui1016]

My cmake may be not wrong, because when I cmake step-6.cc for example, the 
results after cmake is the same as above.

Thanks!

chucui

在 2018年5月27日星期日 UTC+8下午2:41:00，chucu...@gmail.com写道：
>
> Hi, everyone,
>
> I am trying to implement some easier case by using deal.ii, and I have a 
> code and debug it, after cmake and when make it, the results on the 
> Terminal are like this:
>
>
> 
> I type: yum list | grep libgfortran
> the result is:
>
>
> 
> I don't know what happen in my machine, but when I run the examples of 
> deal.ii, every thing is OK.
>
> What should I do next, and how can I debug it successfully? Thank you too 
> much!
>
> chucui
>

-- 
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: when make, something about libgfortran may be wrong

[chucui1016]

After cmake -DDEAL_II_DIR=~/dealii-8.5.1/deal.II, the results is







Maybe my cmake is not successful

Thanks!
chucui


在 2018年5月27日星期日 UTC+8下午2:41:00，chucu...@gmail.com写道：
>
> Hi, everyone,
>
> I am trying to implement some easier case by using deal.ii, and I have a 
> code and debug it, after cmake and when make it, the results on the 
> Terminal are like this:
>
>
> 
> I type: yum list | grep libgfortran
> the result is:
>
>
> 
> I don't know what happen in my machine, but when I run the examples of 
> deal.ii, every thing is OK.
>
> What should I do next, and how can I debug it successfully? Thank you too 
> much!
>
> chucui
>

-- 
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: when make, something about libgfortran may be wrong

[chucui1016]

And this is my code:
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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

#include 
#include 
#include 
#define pi 3.1415926

namespace Case1
{
  using namespace dealii;


  template 
  class PhaseEquation
  {
  public:
PhaseEquation ();
void run ();

  private:
void setup_system ();
void solve_phi_star ();
void assemble_system ();
//void output_results () const;
//double compute_errors () const; //i need to look for this example

Triangulation   triangulation;
FE_Qfe;
DoFHandler  dof_handler;

ConstraintMatrix constraints;

SparsityPattern  sparsity_pattern;
SparseMatrix mass_matrix;
SparseMatrix laplace_matrix;
SparseMatrix matrix_phi;
SparseMatrix mass_bar_matrix;

Vector   solution_phi, solution_phi_star;
Vector   old_solution_phi;
Vector   system_rhs;

double time_step, time;
const double eps;//, s;
unsigned int timestep_number;
// const double eps;
  };


  template 
  PhaseEquation::PhaseEquation () :
fe (1),
dof_handler (triangulation),
time_step (1./64),
time (time_step),
timestep_number (1),
eps (0.03)
  {}


  template 
  class Solution : public Function
 //   protected SolutionBase
  {
  public:
Solution () : Function() {}

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

virtual Tensor<1,dim> gradient (const Point   &p,
const unsigned int  component = 0) 
const;
  };


  // The actual definition of the values and gradients of the exact solution

  template 
  double Solution::value (const Point   &p,
   const unsigned int) const
  {
double return_value = 0;
/*for (unsigned int i=0; in_source_centers; ++i)
  {
const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
return_value += std::exp(-x_minus_xi.norm_square() /
 (this->width * this->width));
  }*/
/*return 0.5 * (1 - std::tanh((p[0] - (3 / (std::sqrt(2) * eps) * 
this->get_time())) / (2 * std::sqrt(2) * eps)))*/
return return_value;
  }



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

/*for (unsigned int i=0; in_source_centers; ++i)
  {
const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];

// For the gradient, note that its direction is along (x-x_i), so we
// add up multiples of this distance vector, where the factor is 
given
// by the exponentials.
return_value += (-2 / (this->width * this->width) *
 std::exp(-x_minus_xi.norm_square() /
  (this->width * this->width)) *
 x_minus_xi);
  }*/

return_value[0] = 0;
return_value[1] = 0;

return return_value;
  }


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

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

  template 
  double InitialValue::value (const Point  &p,
 const unsigned int component) const
  {
(void) component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
const double eps = 0.03;
return 0.5 * (1 - std::tanh((p[0] - (3 / (std::sqrt(2) * eps) * 0)) / 
(2 * std::sqrt(2) * eps)));
  }
  
  
  // Secondly, we have the right hand side forcing term. Boring as we are, 
we
  // choose zero here as well:
  template 
  class RightHandSide : public Function
  {
  public:
RightHandSide () : Function() {}

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



  template 
  double RightHandSide::value (const Point  &/*p*/,
const unsigned int component) const
  {
(void) component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
  }


  // Finally, we have boundary values for $u$ and $v$. They are as described
  // in the introduction, one being the time derivative of the other:
  template 
  class BoundaryValue : public Function
  {
  public:
BoundaryValue () : Function() {}

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

  template 
  double BoundaryValue::value (const Point &p,

[deal.II] when make, something about libgfortran may be wrong

[chucui1016]

Hi, everyone,

I am trying to implement some easier case by using deal.ii, and I have a 
code and debug it, after cmake and when make it, the results on the 
Terminal are like this:


I type: yum list | grep libgfortran
the result is:


I don't know what happen in my machine, but when I run the examples of 
deal.ii, every thing is OK.

What should I do next, and how can I debug it successfully? Thank you too 
much!

chucui

-- 
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] implement weak formulation of biharmonic equation

[chucui1016]

I want to implement grad(delta_u), but shape_grad_component() is function 
in Class FEValuesBase, and it can't be used in delta_u directly, so I want 
to let delta_u be an object in Class FECaluesBase.

在 2018年2月28日星期三 UTC+8上午8:51:11，Wolfgang Bangerth写道：
>
> On 02/27/2018 05:44 PM, chucu...@gmail.com  wrote: 
> > | 
> >  dealii::FEValuesBase<2, dim>::FEValuesBase() ;{} 
> > 
> > | 
>
> The error message shown corresponds to this line. This line does not 
> correspond to anything C++ would recognize. What are you trying to do 
> here? 
>
> 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] implement weak formulation of biharmonic equation

[chucui1016]

ucan/boost/deal.II/examples/practice-4order-39/practice-4order-39.cc:42:
/home/yucan/boost/deal.II/include/deal.II/fe/fe_values.h:2739:3: note: 
dealii::FEValuesBase::FEValuesBase(const 
dealii::FEValuesBase&) [with int dim = 2; int spacedim = 2]
   FEValuesBase (const FEValuesBase &);
   ^
/home/yucan/boost/deal.II/include/deal.II/fe/fe_values.h:2739:3: note:   
candidate expects 1 argument, 0 provided
/home/yucan/boost/deal.II/include/deal.II/fe/fe_values.h:1474:3: note: 
dealii::FEValuesBase::FEValuesBase(unsigned int, unsigned 
int, dealii::UpdateFlags, const dealii::Mapping&, const 
dealii::FiniteElement&) [with int dim = 2; int spacedim = 2]
   FEValuesBase (const unsigned int n_q_points,
   ^
/home/yucan/boost/deal.II/include/deal.II/fe/fe_values.h:1474:3: note:   
candidate expects 5 arguments, 0 provided
In file included from 
/home/yucan/boost/deal.II/examples/practice-4order-39/practice-4order-39.cc:48:0:
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:114:54: 
error: lvalue required as left operand of assignment
 delta_v.shape_value_component(i,k,d) = 
trace(hessian_phi_v);
  ^
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:115:54: 
error: lvalue required as left operand of assignment
 delta_u.shape_value_component(j,k,d) = 
trace(hessian_phi_u);
  ^
make[2]: *** [CMakeFiles/practice4order39.dir/practice-4order-39.cc.o] 
Error 1
make[1]: *** [CMakeFiles/practice4order39.dir/all] Error 2
make: *** [all] Error 2


I wander if my idea is correct and is the problem initialization？
But I've been confused about the right format. Please give me some help! 
Thank you!

chucui1016


在 2018年2月27日星期二 UTC+8上午11:36:22，Timo Heister写道：
>
> > "make run" have error, and I don't understand what it mean and how to 
> do: 
>
> Did you read the error message, namely: 
>
> > You are requesting information from an 
> FEValues/FEFaceValues/FESubfaceValues 
> > object for which this kind of information has not been computed. What 
> > information these objects compute is determined by the update_* flags 
> you 
> > pass to the constructor. Here, the operation you are attempting requires 
> the 
> >  flag to be set, but it was apparently not specified 
> upon 
> > construction. 
>
> ? I think the message is pretty clear on what is wrong. If you 
> struggle with this, you should probably go back and work through the 
> first tutorials. 
>
> -- 
> 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] implement weak formulation of biharmonic equation

[chucui1016]

Dear Timo,

Thank you very much! And I correct what you say, then "make" can go, but 
"make run" have error, and I don't understand what it mean and how to do:
 

yucan@ubuntu-csrc-b224:~/boost/deal.II/examples/practice-4order-39$ make run
[ 50%] Built target practice4order39
[100%] Run with Debug configuration
DEAL::Element: FE_DGQ<2>(3)
DEAL::Step 0
DEAL::Triangulation 16 cells, 2 levels
DEAL::DoFHandler 256 dofs, level dofs 64 256
DEAL::Assemble matrix


An error occurred in line <4236> of file 
 in function
dealii::Tensor<2, spacedim> dealii::FEValuesBase::shape_hessian_component(unsigned int, unsigned int, unsigned 
int) const [with int dim = 2; int spacedim = 2]
The violated condition was: 
this->update_flags & update_hessians
The name and call sequence of the exception was:
ExcAccessToUninitializedField("update_hessians")
Additional Information: 
You are requesting information from an 
FEValues/FEFaceValues/FESubfaceValues object for which this kind of 
information has not been computed. What information these objects compute 
is determined by the update_* flags you pass to the constructor. Here, the 
operation you are attempting requires the  flag to be set, 
but it was apparently not specified upon construction.


Stacktrace:
---
#0  ./practice4order39: dealii::FEValuesBase<2, 
2>::shape_hessian_component(unsigned int, unsigned int, unsigned int) const
#1  ./practice4order39: void 
dealii::LocalIntegrators::Biharmonic::cell_matrix<2>(dealii::FullMatrix&,
 
dealii::FEValuesBase<2, 2> const&, double)
#2  ./practice4order39: 
practice4order39::MatrixIntegrator<2>::cell(dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&) const
#3  ./practice4order39: std::_Function_handler&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&), std::_Bind::*)(dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&) const> 
(dealii::MeshWorker::LocalIntegrator<2, 2, double> const*, 
std::_Placeholder<1>, std::_Placeholder<2>)> >::_M_invoke(std::_Any_data 
const&, dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)
#4  ./practice4order39: std::function&, dealii::MeshWorker::IntegrationInfo<2, 
2>&)>::operator()(dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&) const
#5  ./practice4order39: void 
dealii::MeshWorker::cell_action, dealii::MeshWorker::DoFInfo<2, 2, double>, 2, 2, 
dealii::TriaActiveIterator, false> > 
>(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>&, std::function&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
std::function&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
std::function&, 
dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
dealii::MeshWorker::LoopControl const&)
#6  ./practice4order39: void dealii::WorkStream::run, std::_Placeholder<3>, std::_Placeholder<2>, 
std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, 
dealii::MeshWorker::LoopControl))(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>&, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, 
dealii::MeshWorker::LoopControl const&)>, std::_Bind, 
dealii::MeshWorker::Assembler::MatrixSimple 
>*))(dealii::MeshWorker::DoFInfoBox<2, dealii::MeshWorker::DoFInfo<2, 2, 
double> > const&, 
dealii::MeshWorker::Assembler::MatrixSimple 
>*)>, 
dealii::TriaActiveIterator, false> >, dealii::MeshWorker::IntegrationInfoBox<2, 2>, 
dealii::MeshWorker::DoFInfoBox<2, dealii::MeshWorker::DoFInfo<2, 2, double> 
> 
>(dealii::TriaActiveIterator, false> > const&, 
dealii::identity, false> > >::type const&, std::_Bind, 
std::_Placeholder<3>, std::_Placeholder<2>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, 
dealii::MeshWorker::LoopControl))(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::Me

Re: [deal.II] implement weak formulation of biharmonic equation

[chucui1016]

Dear Timo,

Thank you very much! And I correct what you say, then "make" can go, but 
"make run" have error, and I don't understand what it mean and how to do:
 

yucan@ubuntu-csrc-b224:~/boost/deal.II/examples/practice-4order-39$ make run
[ 50%] Built target practice4order39
[100%] Run with Debug configuration
DEAL::Element: FE_DGQ<2>(3)
DEAL::Step 0
DEAL::Triangulation 16 cells, 2 levels
DEAL::DoFHandler 256 dofs, level dofs 64 256
DEAL::Assemble matrix


An error occurred in line <4236> of file 
 in function
dealii::Tensor<2, spacedim> dealii::FEValuesBase::shape_hessian_component(unsigned int, unsigned int, unsigned 
int) const [with int dim = 2; int spacedim = 2]
The violated condition was: 
this->update_flags & update_hessians
The name and call sequence of the exception was:
ExcAccessToUninitializedField("update_hessians")
Additional Information: 
You are requesting information from an 
FEValues/FEFaceValues/FESubfaceValues object for which this kind of 
information has not been computed. What information these objects compute 
is determined by the update_* flags you pass to the constructor. Here, the 
operation you are attempting requires the  flag to be set, 
but it was apparently not specified upon construction.


Stacktrace:
---
#0  ./practice4order39: dealii::FEValuesBase<2, 
2>::shape_hessian_component(unsigned int, unsigned int, unsigned int) const
#1  ./practice4order39: void 
dealii::LocalIntegrators::Biharmonic::cell_matrix<2>(dealii::FullMatrix&,
 
dealii::FEValuesBase<2, 2> const&, double)
#2  ./practice4order39: 
practice4order39::MatrixIntegrator<2>::cell(dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&) const
#3  ./practice4order39: std::_Function_handler&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&), std::_Bind::*)(dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&) const> 
(dealii::MeshWorker::LocalIntegrator<2, 2, double> const*, 
std::_Placeholder<1>, std::_Placeholder<2>)> >::_M_invoke(std::_Any_data 
const&, dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)
#4  ./practice4order39: std::function&, dealii::MeshWorker::IntegrationInfo<2, 
2>&)>::operator()(dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&) const
#5  ./practice4order39: void 
dealii::MeshWorker::cell_action, dealii::MeshWorker::DoFInfo<2, 2, double>, 2, 2, 
dealii::TriaActiveIterator, false> > 
>(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>&, std::function&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
std::function&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
std::function&, 
dealii::MeshWorker::DoFInfo<2, 2, double>&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>::CellInfo&)> const&, 
dealii::MeshWorker::LoopControl const&)
#6  ./practice4order39: void dealii::WorkStream::run, std::_Placeholder<3>, std::_Placeholder<2>, 
std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, 
dealii::MeshWorker::LoopControl))(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::MeshWorker::IntegrationInfoBox<2, 2>&, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)> const&, 
dealii::MeshWorker::LoopControl const&)>, std::_Bind, 
dealii::MeshWorker::Assembler::MatrixSimple 
>*))(dealii::MeshWorker::DoFInfoBox<2, dealii::MeshWorker::DoFInfo<2, 2, 
double> > const&, 
dealii::MeshWorker::Assembler::MatrixSimple 
>*)>, 
dealii::TriaActiveIterator, false> >, dealii::MeshWorker::IntegrationInfoBox<2, 2>, 
dealii::MeshWorker::DoFInfoBox<2, dealii::MeshWorker::DoFInfo<2, 2, double> 
> 
>(dealii::TriaActiveIterator, false> > const&, 
dealii::identity, false> > >::type const&, std::_Bind, 
std::_Placeholder<3>, std::_Placeholder<2>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, std::function&, dealii::MeshWorker::DoFInfo<2, 
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&, 
dealii::MeshWorker::IntegrationInfo<2, 2>&)>, 
dealii::MeshWorker::LoopControl))(dealii::TriaActiveIterator, false> >, dealii::MeshWorker::DoFInfoBox<2, 
dealii::MeshWorker::DoFInfo<2, 2, double> >&, 
dealii::Me

[deal.II] implement weak formulation of biharmonic equation

[chucui1016]

Dear all,

I want to solve biharmonic equation problem with interior penalty method 
like step-39, but I need to rewrite weak formulation. In step-39, it use 
functions like  LocalIntegrators::Laplace::cell_matrix 

 and 
so on, which define in /include/deal.II/integrators/laplace.h, so I want to 
rewrite a file named biharmonic.h in the directory: 
/include/deal.II/integrators/, but when I compute integration of (delta u * 
delta v) using this code
#ifndef dealii__integrators_biharmonic_h
#define dealii__integrators_biharmonic_h

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

DEAL_II_NAMESPACE_OPEN

namespace LocalIntegrators
{

  namespace Biharmonic
  {

template
void cell_matrix (
  FullMatrix &M,
  const FEValuesBase &fe,
  const double factor = 1.)
{
  const unsigned int n_dofs = fe.dofs_per_cell;
  const unsigned int n_components = fe.get_fe().n_components();
  
  for (unsigned int k=0; k hessian_phi_v = fe.shape_hessian_component(i, 
k, d);
Mii += dx *
   (hessian_phi_v.trace * hessian_phi_v.trace);

  M(i,i) += Mii;

  for (unsigned int j=i+1; j hessian_phi_u = 
fe.shape_hessian_component(j, k, d);
Mij += dx *
   (hessian_phi_u.trace * hessian_phi_v.trace);

  M(i,j) += Mij;
  M(j,i) += Mij;
}
}
}

Then cmake and make, after make, it shows:


yucan@ubuntu-csrc-b224:~/boost/deal.II/examples/practice-4order-39$ make 
[100%] Building CXX object 
CMakeFiles/practice4order39.dir/practice-4order-39.cc.o
In file included from 
/home/yucan/boost/deal.II/examples/practice-4order-39/practice-4order-39.cc:48:0:
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h: In 
function ‘void 
dealii::LocalIntegrators::Biharmonic::cell_matrix(dealii::FullMatrix&, 
const dealii::FEValuesBase&, double)’:
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:54:25: 
error: ‘hessian_phi_v’ was not declared in this scope
(hessian_phi_v.trace * hessian_phi_v.trace);
 ^
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:64:29: 
error: ‘hessian_phi_u’ was not declared in this scope
(hessian_phi_u.trace * hessian_phi_v.trace);
 ^
make[2]: *** [CMakeFiles/practice4order39.dir/practice-4order-39.cc.o] 
Error 1
make[1]: *** [CMakeFiles/practice4order39.dir/all] Error 2
make: *** [all] Error 2


And I also make debug, it shows:


yucan@ubuntu-csrc-b224:~/boost/deal.II/examples/practice-4order-39$ make 
debug
[100%] Switch CMAKE_BUILD_TYPE to Debug
-- Using the deal.II-8.4.2 installation found at /home/yucan/boost/deal.II
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_query_git_information.cmake
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_pickup_tests.cmake
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_initialize_cached_variables.cmake
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_add_test.cmake
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_setup_target.cmake
-- Include macro 
/home/yucan/boost/deal.II/share/deal.II/macros/macro_deal_ii_invoke_autopilot.cmake
-- Configuring done
-- Generating done
-- Build files have been written to: 
/home/yucan/boost/deal.II/examples/practice-4order-39
Scanning dependencies of target practice4order39
[100%] Building CXX object 
CMakeFiles/practice4order39.dir/practice-4order-39.cc.o
In file included from 
/home/yucan/boost/deal.II/examples/practice-4order-39/practice-4order-39.cc:48:0:
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h: In 
function ‘void 
dealii::LocalIntegrators::Biharmonic::cell_matrix(dealii::FullMatrix&, 
const dealii::FEValuesBase&, double)’:
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:54:25: 
error: ‘hessian_phi_v’ was not declared in this scope
(hessian_phi_v.trace * hessian_phi_v.trace);
 ^
/home/yucan/boost/deal.II/include/deal.II/integrators/biharmonic.h:64:29: 
error: ‘hessian_phi_u’ was not declared in this scope
(hessian_phi_u.trace * hessian_phi_v.trace);
 ^
make[6]: *** [CMakeFiles/practice4order39.dir/practice-4order-39.cc.o] 
Error 1
make[5]: *** [CMakeFiles/practice4order39.dir/all] Error 2
make[4]: *** [all] Error 2
make[3]: *** [CMakeFiles/debug] Error 2
make[2]: *** [CMakeFiles/debug.dir/all] Error 2
make[1]: *** [CMakeFiles/debug.dir/rule] Error 2
make: *** [debug] Error 2



I don't know how to correct this question. So ask for help.

Thanks

--