Re: [deal.II] Wrong sign appears in VectorTools::point_gradient function when dealing with hyperball mesh

2018-02-16 Thread Wolfgang Bangerth 

On 02/16/2018 07:34 PM, Moraad Biagooi wrote:


As you suggested, I tried a 2D mesh. The plots shows that the sign difference
happens at the edge of two neighbor symmetric cells answers. Two gradient
functions may calculate on different cells, so it looks like the signs are 
different.


So it looks to me that there was no problem to begin with.


Great, good to know!
Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Wrong sign appears in VectorTools::point_gradient function when dealing with hyperball mesh

2018-02-16 Thread Wolfgang Bangerth 

On 02/15/2018 07:29 AM, Moraad Biagooi wrote:


I have some problems with the signs of VectorTools::point_gradient. 
Here's an example:


I modified Step-4 of tutorial so that the mesh is a 
GridGenerator::hyperball with half of its boundary_id equal to one.


After solving the problem, a function is called that prints 
VectorTools::point_gradient function versus a simple 1st order gradient 
at some points inside the mesh.


The result is almost the same but there's sign difference at some 
points.


Do you know which function is correct? For example, have you taken a 
look at the solution  at that point and determined visually what the 
gradient is, and then compared it with what the two ways of computing 
things return? (This may be easier to do in 2d.)


There is of course always the possibility of a bug, but it would 
surprise me if there was one in the function that computes the gradient 
at a point.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Wrong sign appears in VectorTools::point_gradient function when dealing with hyperball mesh

2018-02-15 Thread Moraad Biagooi 
Dear all,

I have some problems with the signs of VectorTools::point_gradient. Here's 
an example:

I modified Step-4 of tutorial so that the mesh is a 
GridGenerator::hyperball with half of its boundary_id equal to one.

After solving the problem, a function is called that prints  
VectorTools::point_gradient function versus a simple 1st order gradient at 
some points inside the mesh.

The result is almost the same but there's sign difference at some points. 
For example:

correct sign:
 -0.246203 -0.0555845 0.0532539
-0.246227 -0.0555722 0.0532624

wrong sign:
-0.249757 -0.059087 0.0550266
-0.24978 -0.0590716 -0.0550173

correct sign:
 -0.246203 -0.0555845 -0.0532539
-0.246227 -0.0555722 -0.0532454

correct sign:
 -0.242905 -0.0526048 -0.0516154
-0.242929 -0.0525951 -0.0516075

wrong sign:
-0.216536 -0.144069 0.0495166
-0.216534 -0.144071 -0.0495091

Can someone check this out?
I attached the code.

Best regards,
Morad Biagooi

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


using namespace dealii;


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

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void output_results () const;
  void calculate_field () const;

  Triangulation   triangulation;
  FE_Qfe;
  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  Vector   solution;
  Vector   system_rhs;
};



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

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



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

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






template 
double RightHandSide::value (const Point ,
  const unsigned int /*component*/) const
{
  double return_value = 0.0;
  for (unsigned int i=0; i center (0,0,0); 
const double  radius = 1.0;
GridGenerator::hyper_ball 	(triangulation, center,		radius);
 
static const SphericalManifold<3> manifold_description(center);
triangulation.set_manifold (0, manifold_description);

   
triangulation.refine_global(2);

for (typename Triangulation<3,3>::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 (cell->face(f)->center()(0) < 0)
cell->face(f)->set_boundary_id(1);
}

  }
}
  }
  

  std::cout << "   Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< "   Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
}


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

  std::cout << "   Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;

  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern (dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);

  system_matrix.reinit (sparsity_pattern);

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



template 
void Step4::assemble_system ()
{
  QGauss  quadrature_formula(2);


  const RightHandSide right_hand_side;


  FEValues fe_values (fe, quadrature_formula,
   update_values   | update_gradients |
   update_quadrature_points | update_JxW_values);


  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
  const unsigned int   n_q_points= quadrature_formula.size();

  FullMatrix   cell_matrix (dofs_per_cell, dofs_per_cell);
  Vector   cell_rhs (dofs_per_cell);