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 <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <deal.II/base/logstream.h>
using namespace dealii;
template <int dim>
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<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
double return_value = 0.0;
for (unsigned int i=0; i<dim; ++i)
return_value += 4.0 * std::pow(p(i), 4.0);
return return_value;
}
template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
return p.square();
}
template <int dim>
Step4<dim>::Step4 ()
:
fe (1),
dof_handler (triangulation)
{}
template <int dim>
void Step4<dim>::make_grid ()
{
if (dim==3) {
const Point<3> 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<GeometryInfo<3>::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 <int dim>
void Step4<dim>::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 <int dim>
void Step4<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
const RightHandSide<dim> right_hand_side;
FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += (fe_values.shape_grad (i, q_index) *
fe_values.shape_grad (j, q_index) *
fe_values.JxW (q_index));
cell_rhs(i) += (fe_values.shape_value (i, q_index) *
right_hand_side.value (fe_values.quadrature_point (q_index)) *
fe_values.JxW (q_index));
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
BoundaryValues<dim>(),
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim>
void Step4<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> solver (solver_control);
solver.solve (system_matrix, solution, system_rhs,
PreconditionIdentity());
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
}
template <int dim>
void Step4<dim>::output_results () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output (dim == 2 ?
"solution-2d.vtk" :
"solution-3d.vtk");
data_out.write_vtk (output);
}
template <int dim>
void Step4<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
make_grid();
setup_system ();
assemble_system ();
solve ();
output_results ();
calculate_field ();
}
template <int dim>
void Step4<dim>::calculate_field () const {
const double dx = 0.1;
const double dy = 0.1;
const double dz = 0.1;
for (double x = -1; x < 1; x+= dx) {
for (double y = -1; y < 1; y+= dy) {
for (double z = -1; z < 1; z+= dz) {
if ((x*x + y*y + z*z) > .8) continue;
dealii::Point<3> r (x,y,z);
auto f_dr = VectorTools::point_gradient (dof_handler, solution,r);
//------------------------------------
const double inc = 1e-3;
auto rx = r + dealii::Point<3> (inc,0,0);
auto ry = r + dealii::Point<3> (0,inc,0);
auto rz = r + dealii::Point<3> (0,0,inc);
auto f_r0 = VectorTools::point_value (dof_handler, solution,r);
auto f_rx = VectorTools::point_value (dof_handler, solution,rx);
auto f_ry = VectorTools::point_value (dof_handler, solution,ry);
auto f_rz = VectorTools::point_value (dof_handler, solution,rz);
auto f_dx = (f_rx - f_r0)/inc;
auto f_dy = (f_ry - f_r0)/inc;
auto f_dz = (f_rz - f_r0)/inc;
auto f_dr_me = dealii::Point<3> (f_dx, f_dy, f_dz);
//------------------------------------
if (f_dr[0]*f_dr_me[0] < 0 || f_dr[1]*f_dr_me[1] < 0 || f_dr[2]*f_dr_me[2] < 0 )
std::cout << "wrong sign:\n" << f_dr << "\n" << f_dr_me << "\n\n";
else {
std::cout << "correct sign:\n " << f_dr << "\n" << f_dr_me << "\n\n";
}
}
}
}
}
int main ()
{
{
Step4<3> laplace_problem_3d;
laplace_problem_3d.run ();
}
return 0;
}
