Hi all.
I have a question on the solution initialization for iterative methods, in
vector valued problem.
For example, I want to solve 2D flow, where my solution consist of two
block, one is for velocities and the other is for pressure.
First, I construct my initial solution as below
*-Declaration *
...BlockVector<double> previous_solution;
.......
*-Block build*
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_u, n_u);
dsp.block(1,0).reinit (n_p, n_u);
dsp.block(0,1).reinit (n_u, n_p);
dsp.block(1,1).reinit (n_p, n_p);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false
);
sparsity_pattern.copy_from (dsp);
}
*-reinitialization *
previous_solution.reinit (2);
previous_solution.block(0).reinit (n_u);
previous_solution.block(1).reinit (n_p);
previous_solution.collect_sizes ();
After that , I wanted to set my x-componet of velocities by following
command lines.
*VectorTools::interpolate(dof_handler,initialize_solution<dim>(),previous_solution.block(0));*
However this function does not work because my dof_handler has more degree
of freedom. (Degree of freedom of all velocities components + pressure
components)
So I met error message which says
The violated condition was:
vec.size() == dof.n_dofs()
The name and call sequence of the exception was:
ExcDimensionMismatch (vec.size(), dof.n_dofs())
Additional Information:
Dimension 162 not equal to 187
It seems that I have to extract dofs that are only for velocity component
to initialize my u_x velocities.
Is there any way to do this?
Or is there any other way to initialize my particular solution components
only....?
&&FYI initialize_solution<dim> function looks like as follow.
template <int dim>
class initialize_solution : public Function<dim>
{
public:
initialize_solution () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
initialize_solution<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
double x=p[0];
double y=p[1];
double height = 1e-3; // unit of [m]
double return_value;
return_value=std::sin(y*Pi/height);
if (component == 0) // impose only in x direction
return return_value;
else
return 0;
}
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
