how to use post_processor for FE

Jaekwang Kim

unread,

Jul 19, 2017, 6:09:00 PM7/19/17

to deal.II User Group

Hi, all

I was practicing using post_processor, referring step-29 and step 31.

If I want to use post_processor for vector-valued problem, where my solution is (u,v,p)

and I want to compute magnitude of velocity, |u^2+v^2| ignoring pressure parts.

i wrote as step 29;

template <int dim>

class Magnitude : public DataPostprocessorScalar<dim>

{

public:

Magnitude ();

virtual

void

compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,

const std::vector<std::vector<Tensor<1, dim> > > &duh,

const std::vector<std::vector<Tensor<2, dim> > > &dduh,

const std::vector<Point<dim> > &normals,

const std::vector<Point<dim> > &evaluation_points,

std::vector<Vector<double> > &computed_quantities) const;

};

template <int dim>

Magnitude<dim>::Magnitude ()

:

DataPostprocessorScalar<dim> ("ShearRate", update_values)

{}

// Remember that this function may only use data for which the

// respective update flag is specified by

template <int dim>

void

Magnitude<dim>::compute_derived_quantities_vector (

const std::vector<Vector<double> > &uh,

const std::vector<std::vector<Tensor<1, dim> > > & /*duh*/,

const std::vector<std::vector<Tensor<2, dim> > > & /*dduh*/,

const std::vector<Point<dim> > & /*normals*/,

const std::vector<Point<dim> > & /*evaluation_points*/,

std::vector<Vector<double> > &computed_quantities

) const

{

for (unsigned int i=0; i<computed_quantities.size(); i++)

{

Assert(computed_quantities[i].size() == 1,

ExcDimensionMismatch (computed_quantities[i].size(), 1));

Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));

computed_quantities[i](0) = std::sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));

}

and I used it as wrote as step 29;

template <int dim>

void

StokesProblem<dim>::output_results (const unsigned int refinement_cycle) const

{

std::vector<std::string> solution_names (dim, "velocity");

solution_names.push_back ("pressure");

std::vector<DataComponentInterpretation::DataComponentInterpretation>

data_component_interpretation

(dim, DataComponentInterpretation::component_is_part_of_vector);

data_component_interpretation

.push_back (DataComponentInterpretation::component_is_scalar);

DataOut<dim> data_out;

data_out.attach_dof_handler (dof_handler);

data_out.add_data_vector (solution, solution_names,

DataOut<dim>::type_dof_data,

data_component_interpretation);

// added

Magnitude<dim> magnitude;

data_out.add_data_vector (solution.block(0), magnitude); //because I only want to process velocity solution

// added - end ;

data_out.add_data_vector (cellwise_shear_rate,"Shear_rate");

data_out.add_data_vector (cellwise_viscosity,"Viscosity");

data_out.build_patches ();

std::ostringstream filenameeps;

filenameeps << "solution-"<< Utilities::int_to_string (refinement_cycle, 2)<< ".vtk";

std::ofstream output (filenameeps.str().c_str());

data_out.write_vtk (output);

}

but I get error like

An error occurred in line <1022> of file </Users/kimjaekwang/dealii-8.4.1/source/numerics/data_out_dof_data.cc> in function

void dealii::DataOut_DoFData<dealii::DoFHandler<2, 2>, 2, 2>::add_data_vector(const VectorType &, const DataPostprocessor<DoFHandlerType::space_dimension> &) [VectorType = dealii::Vector<double>]

The violated condition was:

vec.size() == dofs->n_dofs()

The name and call sequence of the exception was:

Exceptions::DataOut::ExcInvalidVectorSize (vec.size(), dofs->n_dofs(), dofs->get_triangulation().n_active_cells())

Additional Information:

The vector has size 238 but the DoFHandler object says that there are 274 degrees of freedom and there are 24 active cells. The size of your vector needs to be either equal to the number of degrees of freedom, or equal to the number of active cells.

because my dof_handler has doc for pressure, size does not match...

how can I avoid this problem ?

Wolfgang Bangerth

unread,

Jul 19, 2017, 6:23:36 PM7/19/17

to dea...@googlegroups.com

On 07/19/2017 04:09 PM, Jaekwang Kim wrote:
>
> because my dof_handler has doc for pressure, size does not match...
>
> how can I avoid this problem ?

You just need to pass the *entire* solution vector, not just the first
block. Otherwise, DataOut cannot find out which elements of the vector
correspond to the degrees of freedom enumerated in your DoFHandler.

You will also have to change this:

Magnitude<dim>::compute_derived_quantities_vector ([...])
{
[...]

Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));

[...]

The size will now be 3 (or probably dim+1), even though you are only
using the zeroth and first component of uh in computing the magnitude.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

Jaekwang Kim

unread,

Jul 19, 2017, 7:00:20 PM7/19/17

to deal.II User Group, bang...@colostate.edu

Thanks, !

I tried in a way you advised. so, now i tried.

template <int dim>

void Magnitude<dim>::compute_derived_quantities_vector (

const std::vector<Vector<double> > &uh,

const std::vector<std::vector<Tensor<1, dim> > > & /*duh*/,

const std::vector<std::vector<Tensor<2, dim> > > & /*dduh*/,

const std::vector<Point<dim> > & /*normals*/,

const std::vector<Point<dim> > & /*evaluation_points*/,

std::vector<Vector<double> > &computed_quantities) const

{

for (unsigned int i=0; i<computed_quantities.size(); i++)

{ Assert(computed_quantities[i].size() == 1,

ExcDimensionMismatch (computed_quantities[i].size(), 1));

Assert(uh[i].size() == 3, ExcDimensionMismatch (uh[i].size(), 3));

computed_quantities[i](0) = std::sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));

}

and I use it as ...

template <int dim>

void

StokesProblem<dim>::output_results (const unsigned int refinement_cycle) const

{ [...]

DataOut<dim> data_out;

data_out.attach_dof_handler (dof_handler);

data_out.add_data_vector (solution, solution_names,

DataOut<dim>::type_dof_data,

data_component_interpretation);

// added

Magnitude<dim> magnitude;

data_out.add_data_vector (solution, magnitude); //here solution is block-vector solution

// added - end ;

data_out.build_patches ();

std::ostringstream filenameeps;

filenameeps << "solution-"<< Utilities::int_to_string (refinement_cycle, 2)<< ".vtk";

std::ofstream output (filenameeps.str().c_str());

data_out.write_vtk (output);

[...]

}

but now I receive a error message as...where I don't have clue what's wrong now...

do you have any idea of it ?

An error occurred in line <114> of file </Users/kimjaekwang/dealii-8.4.1/source/base/subscriptor.cc> in function

virtual dealii::Subscriptor::~Subscriptor()

The violated condition was:

counter == 0

The name and call sequence of the exception was:

ExcInUse (counter, object_info->name(), infostring)

Additional Information:

Object of class 9MagnitudeILi2EE is still used by 1 other objects.

(Additional information: <none>)

See the entry in the Frequently Asked Questions of deal.II (linked to from http://www.dealii.org/) for a lot more information on what this error means and how to fix programs in which it happens.

Wolfgang Bangerth

unread,

Jul 19, 2017, 7:02:14 PM7/19/17

to dea...@googlegroups.com

On 07/19/2017 05:00 PM, Jaekwang Kim wrote:
> /See the entry in the Frequently Asked Questions of deal.II (linked to

> from http://www.dealii.org/) for a lot more information on what this

> error means and how to fix programs in which it happens./

Have you looked here? :-)

(The answer is that the `Magnitude` object needs to live at least as
long as the `DataOut` object. You should declare it *before* the
`DataOut` object.)

Cheers

Jaekwang Kim

unread,

Jul 20, 2017, 10:02:47 AM7/20/17

to deal.II User Group, bang...@colostate.edu

Thanks!!

By the way, can I ask another?

Actually, the quantity that I'm going to calculate is \dot(gamma),

More explicitly...I attach one figure

So I modified post processor as,

template <int dim>

void Magnitude<dim>::compute_derived_quantities_vector (

const std::vector<Vector<double> > & uh,

const std::vector<std::vector<Tensor<1, dim> > > & duh,

const std::vector<std::vector<Tensor<2, dim> > > & /*dduh*/,

const std::vector<Point<dim> > & /*normals*/,

const std::vector<Point<dim> > & evaluation_points,

std::vector<Vector<double> > & computed_quantities) const

{

//const unsigned int n_quadrature_points = uh.size();

for (unsigned int i=0; i<computed_quantities.size(); i++)

{ Assert(computed_quantities[i].size() == 1,

ExcDimensionMismatch (computed_quantities[i].size(), 1));

Assert(uh[i].size() == 3, ExcDimensionMismatch (uh[i].size(), 3));

Tensor<2,dim> grad_u;

for (unsigned int d=0; d<dim; ++d)

grad_u[d] = duh[i][d];

std::cout << "grad_u: " << grad_u[0][0] << " , " << grad_u[0][1] << std::end; // to check whether it actually gets value

const SymmetricTensor<2,dim> shear_rate_tensor = symmetrize (grad_u);

computed_quantities[q](0) = shear_rate_tensor * shear_rate_tensor + std::pow( (uh[q](0) /evaluation_points[q](0)),2) );

}

but I got problem in two ways

1. the 'grad_u' vector I actually receiving is all zero vector.... I checked and it always gives me (dudx,dudy)=(0,0)

2. I need to use p[0] value, and I believe 'evaluation_points[q](0)' is what I need.

but when I extract value I am receiving error as...

grad_u: 0 , 0

make[3]: *** [CMakeFiles/run] Segmentation fault: 11

make[2]: *** [CMakeFiles/run.dir/all] Error 2

make[1]: *** [CMakeFiles/run.dir/rule] Error 2

make: *** [run] Error 2

Thanks...

Jaekwang

Wolfgang Bangerth

unread,

Jul 29, 2017, 7:46:28 AM7/29/17

to Jaekwang Kim, deal.II User Group

On 07/20/2017 08:02 AM, Jaekwang Kim wrote:
>
> 1. the 'grad_u' vector I actually receiving is all zero vector.... I checked
> and it always gives me (dudx,dudy)=(0,0)
> 2. I need to use p[0] value, and I believe 'evaluation_points[q](0)' is what I
> need.
> but when I extract value I am receiving error as...

Are you passing the update_gradients and update_quadrature_points flags from
your class to make sure that these two fields are actually filled?

Best
W.

Reply all

Reply to author

Forward

how to use post_processor for FE_System

Jaekwang Kim

Wolfgang Bangerth

Jaekwang Kim

Wolfgang Bangerth

Jaekwang Kim

Wolfgang Bangerth