|
// The vector I want to output is the solution for vector-valued problem.
// I want to output the cell-average value for the 0th component.
Vector<double>modi (solu);
Quadrature<dim>qq(fe.get_unit_support_points());
FEValues<dim>fv(fe,qq,update_values |update_quadrature_points);
std::vector<bool>selected_dofs (dof_handler.n_dofs());
// comp[0] is the 0th component of
std::vector<FEValuesExtractors::Scalar> for the system
DoFTools::extract_dofs(dof_handler,fe.component_mask(comp[0]),selected_dofs);
std::vector<types::global_dof_index>local_to_global_dof (fe.dofs_per_cell);
// modify the 0th component
for(typenameDoFHandler<dim>::active_cell_iterator
cell=dof_handler.begin_active();cell!=dof_handler.end();++cell){
fv.reinit(cell);
std::vector<double>local_sol(qq.size());
fv[comp[0]].get_function_values(solu_out,local_sol);
doubleavg
=std::accumulate(local_sol.begin(),local_sol.end(),0.)/local_sol.size();
cell->get_dof_indices (local_to_global_dof);
for(autoi=0;i<fe.dofs_per_cell;++i){
if(selected_dofs[local_to_global_dof[i]]){
solu_out[local_to_global_dof[i]]=avg;
}
}
}
|
This code gives an approximation in that it takes the average of the
values of the degrees of freedom on each cell, but that is not equal to
the spatial average
1/|K| \int_K u_h(x) dx
that you probably want to compute. That's because if you expand u_h on
cell K into
u_h(x) = \sum_i U_i \phi_i(x)
then the *real* average is given by
\sum_i U_i (1/|K| \int_K phi_i(x) dx)
but in general you do not have that
1/|K| \int_K phi_i(x) dx == const = 1/dofs_per_cell
(This is true for Q1 elements, but not for higher order elements.)
Consequently, you probably want to use a proper integration formula.
Maybe something like this:
QGauss<dim> quadrature(degree + 2);
FEValues<dim> fe_values (...);
Vector<double> cellwise_averages (triangulation.n_active_cells());
for (auto cell : ...)
{
fe_values.reinit (cell);
fe_values.get_function_values (solution, local_values);
double average = 0;
double cell_volume = 0;
for (q=0.....)
{
average += local_values[q] * fe_values.JxW(q);
cell_volume += fe_values.JxW(q);
}
average /= cell_volume;
cellwise_averages(cell->active_cell_index()) = average;
}
You can then attach this vector with as many entries as there are cells
to a DataOut and output it, or do whatever else you need to do with it!
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.
