elasticity/residual.h

Go to the documentation of this file.

/**********************************************************************
 *
 * Copyright Guido Kanschat, 2010, 2012, 2013
 *
 **********************************************************************/

#ifndef __elasticity_residual_h
#define __elasticity_residual_h

#include <amandus/elasticity/integrators.h>
#include <amandus/integrator.h>
#include <deal.II/integrators/divergence.h>
#include <deal.II/integrators/elasticity.h>
#include <deal.II/integrators/grad_div.h>
#include <deal.II/integrators/l2.h>
#include <deal.II/integrators/laplace.h>
#include <elasticity/parameters.h>

using namespace dealii::MeshWorker;

namespace Elasticity
{
template <int dim, int force_type = 0>
class Residual : public AmandusIntegrator<dim>
{
public:
  Residual(const Parameters& par, const dealii::Function<dim>& bdry,
           const std::set<unsigned int>& dirichlet = std::set<unsigned int>());

  virtual void cell(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info) const;
  virtual void boundary(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info) const;
  virtual void face(DoFInfo<dim>& dinfo1, DoFInfo<dim>& dinfo2, IntegrationInfo<dim>& info1,
                    IntegrationInfo<dim>& info2) const;

private:
  dealii::SmartPointer<const Parameters, class Residual<dim>> parameters;
  dealii::SmartPointer<const dealii::Function<dim>, class Residual<dim>> boundary_values;
  std::set<unsigned int> dirichlet_boundaries;
};

//----------------------------------------------------------------------//

template <int dim, int force_type>
Residual<dim, force_type>::Residual(const Parameters& par, const dealii::Function<dim>& bdry,
                                    const std::set<unsigned int>& dirichlet)
  : parameters(&par)
  , boundary_values(&bdry)
  , dirichlet_boundaries(dirichlet)
{
  this->input_vector_names.push_back("Newton iterate");
}

template <int dim, int force_type>
void
Residual<dim, force_type>::cell(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info) const
{
  Assert(info.values.size() >= 1, dealii::ExcDimensionMismatch(info.values.size(), 1));
  Assert(info.gradients.size() >= 1, dealii::ExcDimensionMismatch(info.values.size(), 1));

  const double mu = parameters->mu;
  const double lambda = parameters->lambda;

  if (force_type == 1)
  {
    std::vector<std::vector<double>> force(
      dim, std::vector<double>(info.fe_values(0).n_quadrature_points, 0.));
    std::fill(force[dim - 1].begin(), force[dim - 1].end(), 9.80665);
    dealii::LocalIntegrators::L2::L2(dinfo.vector(0).block(0), info.fe_values(0), force);
  }

  if (parameters->linear)
  {
    dealii::LocalIntegrators::Elasticity::cell_residual(
      dinfo.vector(0).block(0),
      info.fe_values(0),
      dealii::make_slice(info.gradients[0], 0, dim),
      2. * mu);
    dealii::LocalIntegrators::GradDiv::cell_residual(dinfo.vector(0).block(0),
                                                     info.fe_values(0),
                                                     dealii::make_slice(info.gradients[0], 0, dim),
                                                     lambda);
  }
  else
  {
    StVenantKirchhoff::cell_residual(dinfo.vector(0).block(0),
                                     info.fe_values(0),
                                     dealii::make_slice(info.gradients[0], 0, dim),
                                     lambda,
                                     mu);
  }
}

template <int dim, int force_type>
void
Residual<dim, force_type>::boundary(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info) const
{
  std::vector<std::vector<double>> null(
    dim, std::vector<double>(info.fe_values(0).n_quadrature_points, 0.));

  boundary_values->vector_values(info.fe_values(0).get_quadrature_points(), null);

  const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
  if (dirichlet_boundaries.count(dinfo.face->boundary_id()) != 0)
  {
    dealii::LocalIntegrators::Elasticity::nitsche_residual(
      dinfo.vector(0).block(0),
      info.fe_values(0),
      dealii::make_slice(info.values[0], 0, dim),
      dealii::make_slice(info.gradients[0], 0, dim),
      null,
      dealii::LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg),
      2. * parameters->mu);
  }
}

template <int dim, int force_type>
void
Residual<dim, force_type>::face(DoFInfo<dim>& dinfo1, DoFInfo<dim>& dinfo2,
                                IntegrationInfo<dim>& info1, IntegrationInfo<dim>& info2) const
{
  const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
  dealii::LocalIntegrators::Elasticity::ip_residual(
    dinfo1.vector(0).block(0),
    dinfo2.vector(0).block(0),
    info1.fe_values(0),
    info2.fe_values(0),
    dealii::make_slice(info1.values[0], 0, dim),
    dealii::make_slice(info1.gradients[0], 0, dim),
    dealii::make_slice(info2.values[0], 0, dim),
    dealii::make_slice(info2.gradients[0], 0, dim),
    dealii::LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg),
    2. * parameters->mu);
}
}

#endif