maxwell/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_curl_curl_h
#define __matrix_curl_curl_h

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

using namespace dealii;
using namespace LocalIntegrators;

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

  std::vector<double> curl_coefficient;
  std::vector<double> mass_coefficient;
};

template <int dim>
Eigen<dim>::Eigen()
  : curl_coefficient(1, 1.)
  , mass_coefficient(1, 0.)
{
  this->use_boundary = false;
  this->use_face = false;
}

template <int dim>
void
Eigen<dim>::cell(dealii::MeshWorker::DoFInfo<dim>& dinfo,
                 dealii::MeshWorker::IntegrationInfo<dim>& info) const
{
  if (dinfo.n_matrices() > 4)
  {
    AssertDimension(dinfo.n_matrices(), 8);
  }
  else
  {
    AssertDimension(dinfo.n_matrices(), 4);
  }
  const unsigned int id = dinfo.cell->material_id();
  AssertIndexRange(id, curl_coefficient.size());
  AssertIndexRange(id, mass_coefficient.size());
  // const double mu = curl_coefficient[id];
  // const double sigma = mass_coefficient[id];

  Maxwell::curl_curl_matrix(dinfo.matrix(0, false).matrix, info.fe_values(0) /*, mu*/);
  // if (sigma != 0.)
  //   L2::mass_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), sigma);
  Divergence::gradient_matrix(dinfo.matrix(1, false).matrix, info.fe_values(1), info.fe_values(0));
  dinfo.matrix(2, false).matrix.copy_transposed(dinfo.matrix(1, false).matrix);
  if (dinfo.n_matrices() > 4)
  {
    L2::mass_matrix(dinfo.matrix(4, false).matrix, info.fe_values(0));
    //  dinfo.matrix(0,false).matrix.add(-8., dinfo.matrix(4,false).matrix);
  }

  //  L2::mass_matrix(dinfo.matrix(7,false).matrix, info.fe_values(1));
}

template <int dim>
void
Eigen<dim>::boundary(dealii::MeshWorker::DoFInfo<dim>& /*dinfo*/,
                     typename dealii::MeshWorker::IntegrationInfo<dim>& /*info*/) const
{
}

template <int dim>
void
Eigen<dim>::face(dealii::MeshWorker::DoFInfo<dim>& /*dinfo1*/,
                 dealii::MeshWorker::DoFInfo<dim>& /*dinfo2*/,
                 dealii::MeshWorker::IntegrationInfo<dim>& /*info1*/,
                 dealii::MeshWorker::IntegrationInfo<dim>& /*info2*/) const
{
}
}
}

#endif