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]
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 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
* 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, 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/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/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 ()
{
QGauss<dim> quadrature_formula(dof_handler.get_fe().degree+1);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEValues<dim> fe_values_scalar (fe_scalar, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
std::vector<Tensor<1, dim>> fe_function_grad(quadrature_formula.size());
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);
fe_values_scalar.get_function_gradients(test,fe_function_grad);
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));
}
cell_rhs(i) += ((phi_i_u*fe_function_grad[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.add_data_vector (solution, solution_names);
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.add_data_vector (size, "solution");
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());
QGauss<dim> quadrature_formula(dof_handler.get_fe().degree+1);
//Quadrature<dim> quadrature_formula(fe.get_unit_support_points());
std::cout << "test: " << fe.get_unit_support_points() < std::endl;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEValues<dim> fe_values_scalar (fe_scalar, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
std::vector<Tensor<1,dim> > local_velocity_values (quadrature_formula.size());
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
Description: Adobe PDF document
fe3.pdf
Description: Adobe PDF document
#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
data1 = pylab.loadtxt('./output.txt')
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$')
pylab.title('Grad')
pylab.legend(loc=1)
pylab.savefig('plot.pdf')
pylab.show()
