Dear all,
I'm testing MatrixFree framework to build Stokes solver. I'm using Stokes
operator copied from deal.ii tests.
1) I've found out that using Stokes operator initialized for selective use:
div.initialize (system_mf_div_storage, std::vector<unsigned int > (1,1),
std::vector<unsigned int > (1,0)); // div is StokesOperator:
MatrixFreeOperators::Base
for multiplication:
BlockVectorType tb1(1),tb2(1);
tb1.block(0)=system_rhs.block(0);
tb2.block(0)=system_rhs.block(1);
div.vmult(tb2, tb1);
Results in error because vmult relays on local_apply that requires vector
with 2 blocks.
2) I've created divergence operator (B from Stokes). However, vmult checks
requires src and dst to have same size and throws an error. I'm not sure if
pre and post processing of constraints will work on non-square matrix.
Michał Wichrowski
The Divergence operator is defined as follows:
template <int dim,int degree_u, int degree_p, typename Number >
class DivergenceOperator : public MatrixFreeOperators::Base<dim,
LinearAlgebra::distributed::Vector<Number> >
{
public:
typedef typename DoFHandler<dim>::active_cell_iterator CellIterator;
typedef LinearAlgebra::distributed::Vector<Number> VectorType;
DivergenceOperator():
MatrixFreeOperators::Base<dim, VectorType >()
{};
virtual void compute_diagonal()
{}
private:
void
local_apply (const MatrixFree<dim,Number> &data,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int,unsigned int> &cell_range)
const;
virtual void apply_add(VectorType &dst,
const VectorType &src) const;
};
template <int dim,int degree_u, int degree_p, typename Number >
void
DivergenceOperator<dim, degree_u,degree_p,Number>
::apply_add (VectorType &dst,
const VectorType &src) const
{
this->data->cell_loop (&DivergenceOperator::local_apply, this, dst,
src);
}
template <int dim,int degree_u, int degree_p, typename Number >
void
DivergenceOperator<dim, degree_u,degree_p,Number>
::local_apply(const MatrixFree<dim,Number> &data,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int,unsigned int> &cell_range) const
{
typedef VectorizedArray<Number> vector_t;
FEEvaluation<dim,degree_u,degree_u+2,dim,Number> velocity (data, 0);
FEEvaluation<dim,degree_p,degree_u+2,1, Number> pressure (data, 1);
for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
{
velocity.reinit (cell);
velocity.read_dof_values (src);
velocity.evaluate (false,true,false);
pressure.reinit (cell);
for (unsigned int q=0; q<velocity.n_q_points; ++q)
{
vector_t div = -velocity.get_divergence(q);
pressure.submit_value (div, q);
}
pressure.integrate (true,false);
pressure.distribute_local_to_global (dst );
}
}
--
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.
