Dear Forum,
I tried to follow the instructions of the DataPostprocessorTensor class to
output the gradient of a vector valued solution. Therefore I tried to adopt
step-8. However, I am getting the following error message:
An error occurred in line <5640> of file
</opt/dealii_build/source/base/data_out_base.cc> in function
void dealii::DataOutBase::write_vtk(const
std::vector<dealii::DataOutBase::Patch<dim, spacedim> >&, const
std::vector<std::__cxx11::basic_string<char> >&, const
std::vector<std::tuple<unsigned int, unsigned int,
std::__cxx11::basic_string<char, std::char_traits<char>,
std::allocator<char> >,
dealii::DataComponentInterpretation::DataComponentInterpretation> >&, const
dealii::DataOutBase::VtkFlags&, std::ostream&) [with int dim = 2; int
spacedim = 2; std::ostream = std::basic_ostream<char>]
The violated condition was:
std::get<3>(nonscalar_data_range) !=
DataComponentInterpretation::component_is_part_of_tensor
Additional information:
You are trying to use functionality in deal.II that is currently not
implemented. In many cases, this indicates that there simply didn't
appear much of a need for it, or that the author of the original code
did not have the time to implement a particular case. If you hit this
exception, it is therefore worth the time to look into the code to
find out whether you may be able to implement the missing
functionality. If you do, please consider providing a patch to the
deal.II development sources (see the deal.II website on how to
contribute).
I am using deal.ii version 9.4 and am wondering what I did wrong here. Do
you have any suggestions?
Best regards
Konrad
--
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.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/a04ecfe4-3518-4253-930e-addb6f619fd4n%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor.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/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <fstream>
#include <iostream>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem{
public:
ElasticProblem();
void run();
private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results() const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,std::vector<Tensor<1, dim>> & values){
AssertDimension(values.size(), points.size());
Assert(dim >= 2, ExcNotImplemented());
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
values[point_n][0] = 1.0;
else
values[point_n][0] = 0.0;
if (points[point_n].norm_square() < 0.2 * 0.2)
values[point_n][1] = 1.0;
else
values[point_n][1] = 0.0;
}
}
template <int dim>
ElasticProblem<dim>::ElasticProblem()
: dof_handler(triangulation)
, fe(FE_Q<dim>(1), dim)
{}
template <int dim>
void ElasticProblem<dim>::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,
Functions::ZeroFunction<dim>(dim),
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), 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);
}
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_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);
std::vector<double> lambda_values(n_q_points);
std::vector<double> mu_values(n_q_points);
Functions::ConstantFunction<dim> lambda(1.), mu(1.);
std::vector<Tensor<1, dim>> rhs_values(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators()){
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
mu.value_list(fe_values.get_quadrature_points(), mu_values);
right_hand_side(fe_values.get_quadrature_points(), rhs_values);
for (const unsigned int i : fe_values.dof_indices()){
const unsigned int component_i = fe.system_to_component_index(i).first;
for (const unsigned int j : fe_values.dof_indices()){
const unsigned int component_j = fe.system_to_component_index(j).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
{
cell_matrix(i, j) +=
(
(fe_values.shape_grad(i, q_point)[component_i] *
fe_values.shape_grad(j, q_point)[component_j] *
lambda_values[q_point])
+
(fe_values.shape_grad(i, q_point)[component_j] *
fe_values.shape_grad(j, q_point)[component_i] *
mu_values[q_point])
+
((component_i == component_j) ?
(fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) *
mu_values[q_point]) :
0)
) *
fe_values.JxW(q_point);
}
}
}
for (const unsigned int i : fe_values.dof_indices()){
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
cell_rhs(i) += fe_values.shape_value(i, q_point) *
rhs_values[q_point][component_i] *
fe_values.JxW(q_point);
}
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 ElasticProblem<dim>::solve(){
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
}
template <int dim>
class StrainPostprocessor : public DataPostprocessorTensor<dim>{
public:
StrainPostprocessor ()
:
DataPostprocessorTensor<dim> ("strain",
update_gradients)
{}
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> &input_data,
std::vector<Vector<double> > &computed_quantities) const override
{
AssertDimension (input_data.solution_gradients.size(),
computed_quantities.size());
for (unsigned int p=0; p<input_data.solution_gradients.size(); ++p)
{
AssertDimension (computed_quantities[p].size(),
(Tensor<2,dim>::n_independent_components));
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
computed_quantities[p][Tensor<2,dim>::component_to_unrolled_index(TableIndices<2>(d,e))]
= (input_data.solution_gradients[p][d][e]
+
input_data.solution_gradients[p][e][d]) / 2;
}
}
};
template <int dim>
void ElasticProblem<dim>::output_results() const{
StrainPostprocessor<dim> strain;
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(dim, DataComponentInterpretation::component_is_part_of_vector);
data_out.add_data_vector (solution,
std::vector<std::string>(dim,"displacement"),
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.add_data_vector(solution, strain);
data_out.build_patches ();
std::ofstream output("solution.vtk");
data_out.write_vtk (output);
}
template <int dim>
void ElasticProblem<dim>::run(){
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << " Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
setup_system();
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve();
output_results();
}
} // namespace Step8
int main()
{
try
{
Step8::ElasticProblem<2> elastic_problem_2d;
elastic_problem_2d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
