assembly.h

Go to the documentation of this file.

#ifndef ASSEMBLY_H
#define ASSEMBLY_H

#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/fe.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fe_base.h"
#include "libmesh/elem_assembly.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/boundary_info.h"

// rbOOmit includes
#include "libmesh/rb_theta.h"
#include "libmesh/rb_assembly_expansion.h"
#include "libmesh/rb_parametrized_function.h"
#include "libmesh/rb_eim_construction.h"
#include "libmesh/rb_eim_theta.h"

// Bring in bits from the libMesh namespace.
// Just the bits we're using, since this is a header.
using libMesh::boundary_id_type;
using libMesh::DenseSubMatrix;
using libMesh::ElemAssembly;
using libMesh::FEInterface;
using libMesh::FEMContext;
using libMesh::Number;
using libMesh::Point;
using libMesh::RBAssemblyExpansion;
using libMesh::RBConstruction;
using libMesh::RBParameters;
using libMesh::RBParametrizedFunction;
using libMesh::RBTheta;
using libMesh::RBThetaExpansion;
using libMesh::RBEIMAssembly;
using libMesh::RBEIMConstruction;
using libMesh::RBEIMEvaluation;
using libMesh::RBEIMTheta;
using libMesh::Real;
using libMesh::RealGradient;
using libMesh::Elem;
using libMesh::FEBase;
using libMesh::subdomain_id_type;
using libMesh::dof_id_type;

// The function we're approximating with EIM
struct Gxyz : public RBParametrizedFunction
{
  unsigned int get_n_components() const override
  {
    return 3;
  }

  virtual std::vector<Number>
  evaluate(const RBParameters & mu,
           const Point & p,
           dof_id_type /*elem_id*/,
           unsigned int /*qp*/,
           subdomain_id_type /*subdomain_id*/,
           const std::vector<Point> & /*p_perturb*/,
           const std::vector<Real> & /*phi_i_qp*/) override
  {
    Real curvature = mu.get_value("curvature");

    return {1. + curvature*p(0), 1. + curvature*p(0), 1./(1. + curvature*p(0))};
  }
};

struct ThetaA0 : RBTheta
{
  virtual Number evaluate(const RBParameters & mu)
  {
    return mu.get_value("kappa") * mu.get_value("Bi");
  }
};

struct AssemblyA0 : ElemAssembly
{
  virtual void boundary_assembly(FEMContext & c)
  {
    std::vector<boundary_id_type> bc_ids;
    c.get_system().get_mesh().get_boundary_info().boundary_ids (&c.get_elem(), c.side, bc_ids);
    for (const auto & bc_id : bc_ids)
      if (bc_id == 1 || bc_id == 2 || bc_id == 3 || bc_id == 4)
        {
          const unsigned int u_var = 0;

          FEBase * side_fe = nullptr;
          c.get_side_fe(u_var, side_fe);

          const std::vector<Real> & JxW_side = side_fe->get_JxW();

          const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();

          // The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

          // Now we will build the affine operator
          unsigned int n_sidepoints = c.get_side_qrule().n_points();

          for (unsigned int qp=0; qp != n_sidepoints; qp++)
            for (unsigned int i=0; i != n_u_dofs; i++)
              for (unsigned int j=0; j != n_u_dofs; j++)
                c.get_elem_jacobian()(i,j) += JxW_side[qp] * phi_side[j][qp]*phi_side[i][qp];

          break;
        }
  }
};

struct ThetaA1 : RBTheta
{
  virtual Number evaluate(const RBParameters & mu)
  {
    return mu.get_value("kappa") * mu.get_value("Bi") * mu.get_value("curvature");
  }
};

struct AssemblyA1 : ElemAssembly
{
  virtual void boundary_assembly(FEMContext & c)
  {
    std::vector<boundary_id_type> bc_ids;
    c.get_system().get_mesh().get_boundary_info().boundary_ids (&c.get_elem(), c.side, bc_ids);
    for (const auto & bc_id : bc_ids)
      if (bc_id == 1 || bc_id == 3) // y == -0.2, y == 0.2
        {
          const unsigned int u_var = 0;

          FEBase * side_fe = nullptr;
          c.get_side_fe(u_var, side_fe);

          const std::vector<Real> & JxW_side = side_fe->get_JxW();

          const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();

          const std::vector<Point> & xyz = side_fe->get_xyz();

          // The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

          // Now we will build the affine operator
          unsigned int n_sidepoints = c.get_side_qrule().n_points();

          for (unsigned int qp=0; qp != n_sidepoints; qp++)
            {
              Real x_hat = xyz[qp](0);

              for (unsigned int i=0; i != n_u_dofs; i++)
                for (unsigned int j=0; j != n_u_dofs; j++)
                  c.get_elem_jacobian()(i,j) += JxW_side[qp] * x_hat * phi_side[j][qp]*phi_side[i][qp];
            }

          break;
        }
  }
};

struct ThetaA2 : RBTheta {
  virtual Number evaluate(const RBParameters & mu)
  {
    return 0.2*mu.get_value("kappa") * mu.get_value("Bi") * mu.get_value("curvature");
  }
};
struct AssemblyA2 : ElemAssembly
{
  virtual void boundary_assembly(FEMContext & c)
  {
    std::vector<boundary_id_type> bc_ids;
    c.get_system().get_mesh().get_boundary_info().boundary_ids (&c.get_elem(), c.side, bc_ids);
    for (const auto & bc_id : bc_ids)
      if (bc_id == 2 || bc_id == 4) // x == 0.2, x == -0.2
        {
          const unsigned int u_var = 0;

          FEBase * side_fe = nullptr;
          c.get_side_fe(u_var, side_fe);

          const std::vector<Real> & JxW_side = side_fe->get_JxW();

          const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();

          // The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

          // Now we will build the affine operator
          unsigned int n_sidepoints = c.get_side_qrule().n_points();

          if (bc_id == 2)
            {
              for (unsigned int qp=0; qp != n_sidepoints; qp++)
                {
                  for (unsigned int i=0; i != n_u_dofs; i++)
                    for (unsigned int j=0; j != n_u_dofs; j++)
                      c.get_elem_jacobian()(i,j) += JxW_side[qp] * phi_side[j][qp]*phi_side[i][qp];
                }
            }

          if (bc_id == 4)
            {
              for (unsigned int qp=0; qp != n_sidepoints; qp++)
                {
                  for (unsigned int i=0; i != n_u_dofs; i++)
                    for (unsigned int j=0; j != n_u_dofs; j++)
                      c.get_elem_jacobian()(i,j) -= JxW_side[qp] * phi_side[j][qp]*phi_side[i][qp];
                }
            }
        }
  }
};

struct ThetaEIM : RBEIMTheta
{
  ThetaEIM(RBEIMEvaluation & rb_eim_eval_in, unsigned int index_in) :
    RBEIMTheta(rb_eim_eval_in, index_in)
  {}

  virtual Number evaluate(const RBParameters & mu)
  {
    return mu.get_value("kappa") * RBEIMTheta::evaluate(mu);
  }
};

struct AssemblyEIM : RBEIMAssembly
{
  AssemblyEIM(RBEIMConstruction & rb_eim_con_in,
              unsigned int basis_function_index_in) :
    RBEIMAssembly(rb_eim_con_in,
                  basis_function_index_in)
  {}

  virtual void interior_assembly(FEMContext & c)
  {
    // PDE variable numbers
    const unsigned int u_var = 0;

    // EIM variable numbers
    const unsigned int Gx_var = 0;
    const unsigned int Gy_var = 1;
    const unsigned int Gz_var = 2;

    FEBase * elem_fe = nullptr;
    c.get_element_fe(u_var, elem_fe);

    const std::vector<Real> & JxW = elem_fe->get_JxW();

    const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();

    // The number of local degrees of freedom in each variable
    const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

    std::vector<Number> eim_values_Gx;
    evaluate_basis_function(c.get_elem().id(),
                            Gx_var,
                            eim_values_Gx);

    std::vector<Number> eim_values_Gy;
    evaluate_basis_function(c.get_elem().id(),
                            Gy_var,
                            eim_values_Gy);

    std::vector<Number> eim_values_Gz;
    evaluate_basis_function(c.get_elem().id(),
                            Gz_var,
                            eim_values_Gz);

    for (unsigned int qp=0; qp != c.get_element_qrule().n_points(); qp++)
      for (unsigned int i=0; i != n_u_dofs; i++)
        for (unsigned int j=0; j != n_u_dofs; j++)
          c.get_elem_jacobian()(i,j) += JxW[qp] * (eim_values_Gx[qp]*dphi[i][qp](0)*dphi[j][qp](0) +
                                                   eim_values_Gy[qp]*dphi[i][qp](1)*dphi[j][qp](1) +
                                                   eim_values_Gz[qp]*dphi[i][qp](2)*dphi[j][qp](2));
  }
};


struct ThetaF0 : RBTheta
{
  virtual Number evaluate(const RBParameters &) { return 1.; }
};

struct AssemblyF0 : ElemAssembly
{

  virtual void interior_assembly(FEMContext & c)
  {
    const unsigned int u_var = 0;

    FEBase * elem_fe = nullptr;
    c.get_element_fe(u_var, elem_fe);

    const std::vector<Real> & JxW = elem_fe->get_JxW();

    const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();

    // The number of local degrees of freedom in each variable
    const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

    // Now we will build the affine operator
    unsigned int n_qpoints = c.get_element_qrule().n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      for (unsigned int i=0; i != n_u_dofs; i++)
        c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
  }
};

struct ThetaF1 : RBTheta
{
  virtual Number evaluate(const RBParameters & mu)
  {
    return mu.get_value("curvature");
  }
};

struct AssemblyF1 : ElemAssembly
{

  virtual void interior_assembly(FEMContext & c)
  {
    const unsigned int u_var = 0;

    FEBase * elem_fe = nullptr;
    c.get_element_fe(u_var, elem_fe);

    const std::vector<Real> & JxW = elem_fe->get_JxW();

    const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();

    const std::vector<Point> & xyz = elem_fe->get_xyz();

    // The number of local degrees of freedom in each variable
    const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

    // Now we will build the affine operator
    unsigned int n_qpoints = c.get_element_qrule().n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Real x_hat = xyz[qp](0);

        for (unsigned int i=0; i != n_u_dofs; i++)
          c.get_elem_residual()(i) += JxW[qp] * (1.*x_hat*phi[i][qp]);
      }
  }
};

struct Ex6InnerProduct : ElemAssembly
{
  // Assemble the Laplacian operator
  virtual void interior_assembly(FEMContext & c)
  {
    const unsigned int u_var = 0;

    FEBase * elem_fe = nullptr;
    c.get_element_fe(u_var, elem_fe);

    const std::vector<Real> & JxW = elem_fe->get_JxW();

    const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();

    // The number of local degrees of freedom in each variable
    const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

    // Now we will build the affine operator
    unsigned int n_qpoints = c.get_element_qrule().n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      for (unsigned int i=0; i != n_u_dofs; i++)
        for (unsigned int j=0; j != n_u_dofs; j++)
          c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
  }
};

// Define an RBThetaExpansion class for this PDE
// The A terms depend on EIM, so we deal with them later
struct Ex6ThetaExpansion : RBThetaExpansion
{

  Ex6ThetaExpansion()
  {
    attach_A_theta(&theta_a0);
    attach_A_theta(&theta_a1);
    attach_A_theta(&theta_a2);
    attach_F_theta(&theta_f0); // Attach the rhs theta
    attach_F_theta(&theta_f1);
  }

  // The RBTheta member variables
  ThetaA0 theta_a0;
  ThetaA1 theta_a1;
  ThetaA2 theta_a2;
  ThetaF0 theta_f0;
  ThetaF1 theta_f1;
};

// Define an RBAssemblyExpansion class for this PDE
// The A terms depend on EIM, so we deal with them later
struct Ex6AssemblyExpansion : RBAssemblyExpansion
{

  Ex6AssemblyExpansion()
  {
    attach_A_assembly(&assembly_a0);
    attach_A_assembly(&assembly_a1);
    attach_A_assembly(&assembly_a2);
    attach_F_assembly(&assembly_f0); // Attach the rhs assembly
    attach_F_assembly(&assembly_f1);
  }

  // The ElemAssembly objects
  AssemblyA0 assembly_a0;
  AssemblyA1 assembly_a1;
  AssemblyA2 assembly_a2;
  AssemblyF0 assembly_f0;
  AssemblyF1 assembly_f1;
};

#endif