# [deal.II] Re: Vector-valued gradient of solution vector

```Dear Daniel,

To clarify what I am searching for, I would like to get the vector-valued
solution into a format that matches a scalar. Either as a
std::vector<Tensor<1,dim>> or a std::vector<Point<dim>>, the format doesn't
matter as long as it matches a scalar dofs. So I can manipulate the vector
components and for example get the size.```
```
I made a test program that I have attached, that takes the gradient of a
gaussian function. At line 354, I do the convergence from a vector valued
solution to get the size of the vector (as I posted before). This was wrong
and you guys said that I shouldn't assume any particular order of the dofs.
Then I tried to this conversion the correct way, see line 365, but I do not
know how to do this conversion the correct way.

The code I have sent you, produces the correct result for fe1 and fe2, i.e.
when I set fe_scalar (2), fe(FE_Q<dim>(2),dim) at line 100. But failed for
fe3. I have attached the result I got for fe2 and fe3 as well as the plot
script I have used to generate these plots.

I would be really grateful if you could help me figure out how to solve
this.

Best,
Joel

On Monday, September 19, 2016 at 11:56:51 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
>> [...]
>>      for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
>>        for (unsigned int i=0; i<dofs_per_cell; ++i)
>>          {
>>         velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>>          }
>>
>> Here you assign the same Tensor<1,dim> at all DoFs multiple times and you
> end with velocity[*]=local_velocity_values[n_q_points]-1.
> I don't think that this is what you want to do, is it? If I understand you
> correctly, you try to assign the magnitude of the velocity in each degree
> of freedom
> to a different scalar Vector that represents a scalar FiniteElement Vector.
> In this case, you need to use a Quadrature object that uses the
> unit_support_points of the scalar FE [1] and your assignment would read
>
> for (unsigned int i=0; i<dofs_per_cell; ++i)
>   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>
> where velocity is a Vector that is related to fe_scalar and
> local_velocity_values has the type  std::vector<Vector<double> >.
>
> Best,
> Daniel
>
> [1]
>

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
---
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
```
```/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2013 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
* 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, 1999
*/

#include <deal.II/grid/tria.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.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/compressed_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparse_direct.h>

#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <math.h>

#include <deal.II/base/logstream.h>

using namespace dealii;

const double pi=3.14159265359;
unsigned int center_id=0;

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;

Triangulation<dim>   triangulation;
FE_Q<dim>            fe_scalar;
FESystem<dim>		   fe;

DoFHandler<dim>      dof_handler_scalar;
DoFHandler<dim>      dof_handler;

SparsityPattern      sparsity_pattern;
SparseMatrix<double> system_matrix;

ConstraintMatrix     constraints;

Vector<double>       solution;
Vector<double>       system_rhs;
Vector<double>       test;
Vector<double>       size;

std::vector<Point<dim>> nodes;
};

template <int dim>
Step4<dim>::Step4 ()
:
fe_scalar (3),
fe (FE_Q<dim>(3),dim),
dof_handler (triangulation),
dof_handler_scalar (triangulation)
{}

template <int dim>
void Step4<dim>::make_grid ()
{
GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
triangulation.refine_global (2);

for (unsigned int step=0; step<2; ++step)
{
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int v=0;v < GeometryInfo<dim>::vertices_per_cell;++v)
{
Point<dim> center;
for (unsigned int i=0; i<dim; ++i)
center(i)=0;

const double distance_from_center = center.distance (cell->vertex(v));
if (std::fabs(distance_from_center) < 1.0)
{
cell->set_refine_flag ();
break;
}
}
}
triangulation.execute_coarsening_and_refinement ();
}

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

dof_handler_scalar.distribute_dofs (fe_scalar);

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

std::cout << "   Number of degrees of freedom for scalar: "
<< dof_handler_scalar.n_dofs()
<< std::endl;

constraints.clear ();

DoFTools::make_hanging_node_constraints (dof_handler,constraints);

constraints.close ();

CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, c_sparsity,constraints,/*keep_constrained_dofs = */false);
sparsity_pattern.copy_from(c_sparsity);

system_matrix.reinit (sparsity_pattern);

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

test.reinit (dof_handler_scalar.n_dofs());

size.reinit (dof_handler_scalar.n_dofs());

nodes.resize(dof_handler_scalar.n_dofs());

MappingQ1<dim> mapping;
DoFTools::map_dofs_to_support_points(mapping,dof_handler_scalar,nodes);

for (unsigned int i=0;i<test.size();i++)
{
Point<dim> center;
for (unsigned int j=0; j<dim; ++j)
center(j)=0;

if (nodes[i].distance(center) < 1e-8)
{
center_id = i;
std::cout << "center id: " << center_id << std::endl;
}
}

std::cout << "fe_degree: " << dof_handler.get_fe().degree << std::endl;
std::cout << "quad_degree: " << dof_handler.get_fe().degree+1 << std::endl;

std::cout << "nodes: " << nodes[center_id] << std::endl;
}

template <int dim>
void Step4<dim>::assemble_system ()
{

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_scalar = dof_handler_scalar.begin_active(),
endc_scalar = dof_handler_scalar.end();

typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();

for (; cell!=endc; ++cell, ++cell_scalar)
{
fe_values.reinit (cell);
fe_values_scalar.reinit (cell_scalar);
cell_matrix = 0;
cell_rhs = 0;

const FEValuesExtractors::Vector velocities (0);

for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> phi_i_u     = fe_values[velocities].value (i, q_index);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> phi_j_u     = fe_values[velocities].value (j, q_index);

cell_matrix(i,j) += ((phi_i_u*phi_j_u)*
fe_values.JxW (q_index));
}
fe_values.JxW (q_index));
}

cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}

template <int dim>
void Step4<dim>::solve ()
{
SolverControl           solver_control (1000, 1e-6);
SolverCG<>              solver (solver_control);

PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);

solver.solve (system_matrix, solution, system_rhs, preconditioner);

std::cout << "   " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;

constraints.distribute (solution);
}

template <int dim>
void Step4<dim>::output_results () const
{
// solution
DataOut<dim> data_out;

// vector
std::vector<std::string> solution_names;
solution_names.push_back ("x");
solution_names.push_back ("y");
solution_names.push_back ("z");

data_out.attach_dof_handler (dof_handler);

data_out.build_patches ();

std::ofstream output ("./solution.vtk");
data_out.write_vtk (output);

//test
DataOut<dim> data_out2;

data_out2.attach_dof_handler (dof_handler_scalar);

data_out2.build_patches ();

std::ofstream output2 ("./size.vtk");
data_out2.write_vtk (output2);

//text file
std::ofstream myfile;
myfile.open ("./output.txt", std::ios::out | std::ios::trunc);

for (unsigned int i=0;i<dof_handler_scalar.n_dofs();i++)
{
if (fabs(nodes[i](1)) < 1e-8 and fabs(nodes[i](2)) < 1e-8)
{
myfile << nodes[i] << "\t" << size[i] << std::endl;
}
}

myfile.close();
}

template <int dim>
void Step4<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;

make_grid();
setup_system ();

Point<dim> center;
for (unsigned int i=0; i<dim; ++i)
center(i)=0;

double factor = 1.0;

for (unsigned int i=0;i<test.size();i++)
{
//Gaussian test function
test(i)= exp(-pow((nodes[i].distance(center)),2.0));
}

assemble_system ();
solve ();

for (unsigned int i=0;i<solution.size();i=i+3)
{
unsigned int x=i;
unsigned int y=i+1;
unsigned int z=i+2;
unsigned int j = i / 3.0;

size(j)=pow(solution(x),2.0)+pow(solution(y),2.0)+pow(solution(z),2.0);
}

/*
std::vector<Tensor<1,dim>>       velocity;
velocity.resize (dof_handler_scalar.n_dofs());

std::cout << "test: " << fe.get_unit_support_points() < std::endl;

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

std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

const FEValuesExtractors::Vector velocities (0);

typename DoFHandler<dim>::active_cell_iterator
cell_scalar = dof_handler_scalar.begin_active(),
endc_scalar = dof_handler_scalar.end();

typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();

for (; cell!=endc; ++cell, ++cell_scalar)
{
fe_values.reinit (cell);
fe_values_scalar.reinit (cell_scalar);
fe_values[velocities].get_function_values (solution,local_velocity_values);

cell_scalar->get_dof_indices (local_dof_indices);

for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
velocity[local_dof_indices[i]] = local_velocity_values[q_index];
}
}

for (unsigned int i = 0;i <size.size()  ;i++ )
{
size(i) = velocity[i]*velocity[i];
}
*/
output_results ();
}

int main ()
{
deallog.depth_console (0);

Step4<3> test;
test.run();

return 0;
}
```

fe2.pdf

fe3.pdf

```#import matplotlib as mpl
#mpl.use('pdf')
import pylab
import numpy as np

pylab.rcParams['axes.linewidth'] = 2
pylab.rcParams['xtick.major.size'] = 5
pylab.rcParams['xtick.major.width'] = 1.5
pylab.rcParams['xtick.minor.size'] = 2.5
pylab.rcParams['xtick.minor.width'] = 1.5
pylab.rcParams['ytick.major.size'] = 5
pylab.rcParams['ytick.major.width'] = 1.5
pylab.rcParams['ytick.minor.size'] = 2.5
pylab.rcParams['ytick.minor.width'] = 1.5

xvals1 = data1[:,0]
yvals1 = data1[:,3]

order1 = np.argsort(xvals1)

xs1 = np.array(xvals1)[order1]
ys1 = np.array(yvals1)[order1]

# analytical soultion
r = np.arange(-2.5, 2.5, 0.01)
s = 4*pow(r,2)*np.exp(-2*pow(r,2))

pylab.plot(xs1,ys1,'r.-',label='dealii',linewidth=2,markersize=10,markeredgewidth=2.0,markerfacecolor='none')

pylab.plot(r,s,'k-',label='exact',linewidth=2,markersize=10,markeredgewidth=2.0,markerfacecolor='none')

pylab.xlabel('X')
pylab.ylabel(r'\$|\nabla(e^{-r^2})|^2\$')