rb_classes.h

Go to the documentation of this file.

#ifndef RB_CLASSES_H
#define RB_CLASSES_H

// local includes
#include "assembly.h"

// rbOOmit includes
#include "libmesh/rb_construction.h"
#include "libmesh/rb_evaluation.h"

// libMesh includes
#include "libmesh/fe_base.h"
#include "libmesh/dof_map.h"

using namespace libMesh;

#ifdef LIBMESH_ENABLE_DIRICHLET

class ElasticityRBEvaluation : public RBEvaluation
{
public:

  ElasticityRBEvaluation(const Parallel::Communicator & comm)
    : RBEvaluation(comm)
  {
    set_rb_theta_expansion(elasticity_theta_expansion);
  }

  virtual Real get_stability_lower_bound() { return 1.; }

  ElasticityThetaExpansion elasticity_theta_expansion;
};


class ElasticityRBConstruction : public RBConstruction
{
public:

  ElasticityRBConstruction (EquationSystems & es,
                            const std::string & name_in,
                            const unsigned int number_in) :
    Parent(es, name_in, number_in),
    u_var(0), v_var(0), w_var(0), var_order(FIRST),
    elasticity_assembly_expansion(*this),
    ip_assembly(*this)
  {}

  virtual ~ElasticityRBConstruction () = default;

  typedef ElasticityRBConstruction sys_type;

  typedef RBConstruction Parent;

  virtual void init_data()
  {
    unsigned int n_components = this->get_mesh().mesh_dimension();

    libmesh_assert_less_equal(n_components, 3);

    u_var = this->add_variable("u", var_order);
    if (n_components > 1)
      v_var = this->add_variable("v", var_order);
    if (n_components > 2)
      w_var = this->add_variable("w", var_order);

    // Generate a DirichletBoundary object
    dirichlet_bc = build_zero_dirichlet_boundary_object();

    // Set the Dirichlet boundary condition
    dirichlet_bc->b.insert(BOUNDARY_ID_MIN_X); // Dirichlet boundary at x=0
    dirichlet_bc->variables.push_back(u_var);
    if (n_components > 1)
      dirichlet_bc->variables.push_back(v_var);
    if (n_components > 2)
      dirichlet_bc->variables.push_back(w_var);

    // Attach dirichlet_bc (must do this _before_ Parent::init_data)
    get_dof_map().add_dirichlet_boundary(*dirichlet_bc);

    Parent::init_data();

    // Set the rb_assembly_expansion for this Construction object
    set_rb_assembly_expansion(elasticity_assembly_expansion);

    // We need to define an inner product matrix for this problem
    set_inner_product_assembly(ip_assembly);
  }

  virtual void init_context(FEMContext & c)
  {
    // For efficiency, we should prerequest all
    // the data we will need to build the
    // linear system before doing an element loop.
    FEBase * elem_fe = nullptr;
    c.get_element_fe(u_var, elem_fe);

    elem_fe->get_JxW();
    elem_fe->get_phi();
    elem_fe->get_dphi();

    FEBase * side_fe = nullptr;
    c.get_side_fe(u_var, side_fe);
    side_fe->get_JxW();
    side_fe->get_phi();
  }

  unsigned int u_var;
  unsigned int v_var;
  unsigned int w_var;
  Order var_order;

  ElasticityAssemblyExpansion elasticity_assembly_expansion;

  InnerProductAssembly ip_assembly;

  std::unique_ptr<DirichletBoundary> dirichlet_bc;
};

#endif // LIBMESH_ENABLE_DIRICHLET

#endif