stokes/eigen.h

Go to the documentation of this file.

/**********************************************************************
 *  Copyright (C) 2011 - 2014 by the authors
 *  Distributed under the MIT License
 *
 * See the files AUTHORS and LICENSE in the project root directory
 **********************************************************************/
#ifndef __matrix_stokes_h
#define __matrix_stokes_h

#include <amandus/integrator.h>
#include <deal.II/integrators/divergence.h>
#include <deal.II/integrators/l2.h>
#include <deal.II/integrators/laplace.h>
#include <deal.II/meshworker/integration_info.h>

using namespace dealii;
using namespace LocalIntegrators;

namespace StokesIntegrators
{
template <int dim>
class Eigen : public AmandusIntegrator<dim>
{
public:
  virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo,
                        MeshWorker::IntegrationInfo<dim>& info) const;
  virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
                    MeshWorker::IntegrationInfo<dim>& info1,
                    MeshWorker::IntegrationInfo<dim>& info2) const;
};

template <int dim>
void
Eigen<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
{
  if (dinfo.n_matrices() > 4)
  {
    AssertDimension(dinfo.n_matrices(), 8);
  }
  else
  {
    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);

  // Eigenvalue mass matrices
  if (dinfo.n_matrices() > 4)
  {
    L2::mass_matrix(dinfo.matrix(4, false).matrix, info.fe_values(0));
    //  L2::mass_matrix(dinfo.matrix(7,false).matrix, info.fe_values(1));
  }
}

template <int dim>
void
Eigen<dim>::boundary(MeshWorker::DoFInfo<dim>& dinfo,
                     typename MeshWorker::IntegrationInfo<dim>& info) const
{
  const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
  Laplace::nitsche_matrix(dinfo.matrix(0, false).matrix,
                          info.fe_values(0),
                          Laplace::compute_penalty(dinfo, dinfo, deg, deg));
}

template <int dim>
void
Eigen<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));
}
}

#endif