MeshWorker and vector valued problem

[Michał Wichrowski]

I dealing with system of equations and I'll be using DG method to solve them.
I'm trying to use MeshWorker for system assembly, but I have no idea how do do it.
I've alredy searched for it in tutorials and MeshWorker namespace documentation.
Michał Wichrowski

[Guido Kanschat]

Dear Michal,

in the simplest case, that is, if you assemble matrices using extractors, there is no difference at all. What I typically do is assembling cell matrices as blocks, in which case you should switch on local renumbering.

Let me attach a code snippet for matrix assembling. First, the global loop. Note the two statements involving dof_handler.block_info(). These trigger the local renumbering.

template <int dim>

void

Application<dim>::assemble_matrix()

{

  matrix = 0.;

  MeshWorker::IntegrationInfoBox<dim> info_box;

  UpdateFlags update_flags = update_values | update_gradients | update_hessians | update_quadrature_points;

  info_box.add_update_flags_all(update_flags);

  info_box.initialize(fe, mapping, &dof_handler.block_info());

  MeshWorker::DoFInfo<dim> dof_info(dof_handler.block_info());

  MeshWorker::Assembler::MatrixSimple<SparseMatrix<double> > assembler;

  assembler.initialize_local_blocks(dof_handler.block_info().local());

  assembler.initialize(matrix);

  assembler.initialize(constraints);

  MeshWorker::integration_loop<dim, dim>(

    dof_handler.begin_active(), dof_handler.end(),

    dof_info, info_box, matrix_integrator, assembler);

}

=========

The matrix_integrator object is derived from MeshWorker::LocalIntegrator (a recent class simplifying the call to the loop). As an example, here are the cell and face function for the Stokes problem:

template <int dim>

void StokesIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const

{

  AssertDimension (dinfo.n_matrices(),4);

  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0));

  Divergence::cell_matrix(dinfo.matrix(2,false).matrix, info.fe_values(0), info.fe_values(1));

  dinfo.matrix(1,false).matrix.copy_transposed(dinfo.matrix(2,false).matrix);

}

template <int dim>

void StokesIntegrator<dim>::face(

  MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,

  MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const

{

  const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();

  Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 

    dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,

    info1.fe_values(0), info2.fe_values(0),

    Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));

}

The numbering of the matrices is row wise in the scheme

A G

D 0

All this uses a single matrix and a single vector, but they can be of block type. I hope it helps.

Best,

Guido

[Michał Wichrowski]

Thanks for answer.
For verifying that I understand you correctly:
I'm modifying step-12 program and adding second equation ( Laplace equation, completely unconnected with first equation).
So my test case will look like
div(beta u) =0
div(grad(p) )=0

after :
-----------
    MeshWorker::Assembler::SystemSimple<SparseMatrix<double>,Vector<double> > assembler;
-----------
I should add:
----------
  assembler.initialize_local_blocks(dof_handler.block_info().local());
----------

And integrate_cell_term function should look like:
====================SOURCE CODE================
template<int dim>
void AdvectionProblem<dim>::integrate_cell_term(DoFInfo& dinfo,
                                     CellInfo& info)
{
   /* LINE CHANGED */
    const FEValuesBase<dim>& fe_v = info.fe_values(0); 

    FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
    const std::vector<double>& JxW =fe_v.get_JxW_values();

    for(unsigned int point =0; point<fe_v.dofs_per_cell; ++point )
    {
        Point<dim> beta;
        beta(0) = -fe_v.quadrature_point(point)(1);
        beta(1) =  fe_v.quadrature_point(point)(0);
        beta /=beta.norm();
        for(unsigned int i =0; i<fe_v.dofs_per_cell; ++i)
        {
            for(unsigned int j = 0; j< fe_v.dofs_per_cell; ++j)
            {
                local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)
                                    *fe_v.shape_value(j,point)
                                    *JxW[point];
            }
        }
    }

    /*LINES ADDED */
    const FEValuesBase<dim> fe_v_2nd = info.fe_values(1);
    FullMatrix<double>& local_matrix_2nd = dinfo.matrix(2).matrix;
   
    for(unsigned int point =0; point<fe_v.dofs_per_cell; ++point )
    {
        for(unsigned int i =0; i<fe_v.dofs_per_cell; ++i)
        {
            for(unsigned int j = 0; j< fe_v.dofs_per_cell; ++j)
            {
                local_matrix(i,j) += fe_v_2nd.shape_grad(i,point)
                                    *fe_v_2nd.shape_grad(j,point)
                                    *JxW[point];
            }
        }
    }
}
==========================================

And I should make similar changes in face and boundary integrators.


What is the numbering of matrixes?
Is it like that:
 0  1  4
 2  3  5 
 6  7  8
(when adding new variable to system previous numbering does not change)
or:
0   1  2
3   4  5
6   7  8

Is it possible to use FEValuesExtractors and assemble matrix like in step-20 program ?
Michał Wichrowski

[Guido Kanschat]

Dear Michał,

It looks good to me, except for the matrix number, which should be 3, not 2, if I understand you correctly.

> 0   1  2
> 3   4  5
> 6   7  8

This is the numbering scheme.

>
> Is it possible to use FEValuesExtractors and assemble matrix like in step-20 program ?

Yes. Just do not initialize with block info.

Then, you have a single matrix like in the onedimensional case.

[Michał Wichrowski]

Dear Guido,
i've managed to make use of FEValuesExtractors. It works ! Thanks for help.
what does following line :

    FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;

?
What does 0 mean?

Thanks for help
Michał Wichrowski

[Guido Kanschat]

Michał,

In your case the zero only indicates that there might be more.

Another model of assembling is that you separately fill several blocks or matrices. In that case, you  would have several matrices.

Best,
Guido

--
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
 
 

[Ernesto Ismail]

Hi,

I'm trying to follow this post to learn how to solve vector valued problems with Meshworker.

Two things, both related to assembler.initialize_local_blocks(dof_handler.block_info().local()); as Mishal implemented in modifying step-12 :

Firstly: This function does not seem to be defined for SystemSimple type objects in Meshworker. Am I missing something in trying to follow Mishal's steps? I note that it is defined for the Matrix and Vector objects however.

Secondly: The comments note that this function has been depreciated. What is the preferred manner to achieve block structures in the manner discussed above?

Lastly: is this still the preferred manner to create systems to solve vector valued problems?

Thanks,

Ernesto