small_strain_mech.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2015:  G-CSC, Goethe University Frankfurt
 * Authors: Raphael Prohl, Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef SMALL_STRAIN_MECH_H_
#define SMALL_STRAIN_MECH_H_
 
#include <stdio.h>
#include <string>
 
#include "bridge/bridge.h"
 
#include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h"
#include "lib_disc/spatial_disc/disc_util/fe_geom.h"
#include "lib_disc/local_finite_element/lagrange/lagrange.h"
#include "lib_disc/quadrature/gauss/gauss_quad.h"
#include "lib_disc/spatial_disc/user_data/data_export.h"
#include "lib_disc/spatial_disc/user_data/data_import.h"
 
#include "material_laws/mat_law_interface.h"
#include "output_writer/mech_output_writer.h"
 
namespace ug{
namespace SmallStrainMechanics{
 
 
template <typename TDomain>
class SmallStrainMechanicsElemDisc
  : public IElemDisc<TDomain>
{
  private:
    typedef IElemDisc<TDomain> base_type;
 
    typedef SmallStrainMechanicsElemDisc<TDomain> this_type;
 
    typedef typename domain_traits<TDomain::dim>::grid_base_object TBaseElem;
 
  public:
    typedef typename base_type::domain_type domain_type;
 
    static const int dim = base_type::dim;
 
    typedef typename base_type::position_type position_type;
 
  public:
    SmallStrainMechanicsElemDisc(const char* functions, const char* subsets);
 
    virtual ~SmallStrainMechanicsElemDisc();
 
    void set_material_law(SmartPtr<IMaterialLaw<TDomain> > spMatLaw)
    { m_spMatLaw = spMatLaw;}
 
    SmartPtr<IMaterialLaw<TDomain> > get_material_law()
    { return m_spMatLaw;}
 
    void set_output_writer(SmartPtr<MechOutputWriter<TDomain> > spOutWriter){
      m_spOutWriter = spOutWriter; m_bOutWriter = true;
      m_spOutWriter->preprocess();
    }
 
    void set_volume_forces(SmartPtr<CplUserData<MathVector<dim>, dim> > user);
    void set_volume_forces(number vel);
    void set_volume_forces(number vel_x, number vel_y);
    void set_volume_forces(number vel_x, number vel_y, number vel_z);
  #ifdef UG_FOR_LUA
      void set_volume_forces(const char* fctName);
  #endif
 
    void set_div_factor(SmartPtr<CplUserData<number, dim> > user);
    void set_div_factor(number val);
#ifdef UG_FOR_LUA
    void set_div_factor(const char* fctName);
#endif
 
    void set_compress_factor(number val);
 
    void set_viscous_forces(SmartPtr<CplUserData<MathVector<dim>, dim> > user0,
                SmartPtr<CplUserData<MathVector<dim>, dim> > user1);
 
    void set_quad_order(const size_t order) {m_quadOrder = order; m_bQuadOrderUserDef = true;}
    int get_quad_order() {return m_quadOrder;}
 
    void set_mass_scale(double val) { m_massScale = val; }
 
    void init_state_variables(const size_t order);
 
    template<typename TElem>
    void contact_forces_elem_ips_avg(LocalVector& locForce, GridObject* side, TElem* elem,
        const MathVector<dim> sideCoPos[], int numElemCorners, const LocalVector& locU,
        std::vector<DoFIndex> vActiveSetLoc);
 
    template <typename TSide, typename TElem>
    void contact_forces_elem_midpoint(LocalVector& locForce, TSide* side, TElem* elem,
        const MathVector<dim> sideCoPos[], const LocalVector& locU,
        std::vector<DoFIndex> vActiveSetLoc);
 
 
    virtual void prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
    {
    //  check number
      if(vLfeID.size() != (size_t)dim)
        UG_THROW("SmallStrainMechanics: Needs exactly "<<dim<<" functions.");
 
      //  check & set order
      for(size_t i = 0; i < vLfeID.size(); ++i)
        if(vLfeID[i].order() < 1)
          UG_THROW("SmallStrainMechanics: Adaptive order not implemented.");
 
    //  remember lfeID;
      m_lfeID = vLfeID[0];
 
    //  set order
      m_order = vLfeID[0].order();
 
    //  update assemble functions
      set_assemble_funcs();
    }
 
 
    std::string config_string() const
    {
      std::stringstream ss;
      ss << "SmallStrainMechanics " << dim << "d ";
      if(dim == 2)
        ss << " [Plain Strain / Ebener Verzerrungszustand]";
      ss << "\n";
      ss << " order of disc scheme = " << m_order
          << ", shape function set " << m_lfeID
          << ", Mass Scale = " << m_massScale
          << "\n";
      if(m_bQuadOrderUserDef)
        ss << " User Defined Quad Order = " << m_quadOrder << "\n";
 
      ss << " Material Configuration: " << ConfigShift(m_spMatLaw->m_materialConfiguration) << "\n";
 
      return ss.str();
    }
 
  private:
    size_t num_fct() const {return dim;}
 
    void set_assemble_funcs();
 
    void register_all_fe_funcs(int order, int quadOrder);
    template <typename TElem, typename TFEGeom>
    void register_fe_func();
 
    template<typename TElem, typename TFEGeom>
    void prep_timestep_elem(const number time, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void prep_elem_loop(const ReferenceObjectID roid, const int si);
 
    template<typename TElem, typename TFEGeom>
    void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void fsh_elem_loop();
 
    template<typename TElem, typename TFEGeom>
    void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    template<typename TElem, typename TFEGeom>
    void fsh_timestep_elem(const number time, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
 
    void update_geo_elem(TBaseElem* elem, DimFEGeometry<dim>& geo);
 
    void print_mat_constants(const number lambda,
        const number mu, const number E, const number v);
 
  protected:
    template <typename TElem, typename TFEGeom>
    void lin_def_pressure(const LocalVector& u,
                std::vector<std::vector<number> > vvvLinDef[],
                const size_t nip);
 
    template <typename TElem, typename TFEGeom>
    void lin_def_volume_forces(const LocalVector& u,
                  std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
                  const size_t nip);
 
 
    template <typename TElem, typename TFEGeom>
    void lin_def_viscous_forces0(const LocalVector& u,
                      std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
                      const size_t nip);
 
    template <typename TElem, typename TFEGeom>
    void lin_def_viscous_forces1(const LocalVector& u,
                          std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
                          const size_t nip);
 
    template <typename TElem, typename TFEGeom>
    void ex_displacement(const LocalVector& u,
               const MathVector<dim> vGlobIP[],
               const MathVector<TFEGeom::Type::dim> vLocIP[],
               const size_t nip,
               MathVector<dim> vValue[],
               bool bDeriv,
               std::vector<std::vector<MathVector<dim> > > vvvDeriv[]);
 
  private:
    SmartPtr<IMaterialLaw<TDomain> > m_spMatLaw;
 
    SmartPtr<MathTensor4<dim,dim,dim,dim> > m_spElastTensor;
 
    SmartPtr<MechOutputWriter<TDomain> > m_spOutWriter;
    bool m_bOutWriter;
    bool m_bMatLawPassedToOutWriter;
 
    int m_order;
 
    LFEID m_lfeID;
 
    bool m_bQuadOrderUserDef;
    int m_quadOrder;
 
    double m_massScale;
 
    DataImport<MathVector<dim>, dim > m_imVolForce;
 
    DataImport<number, dim> m_imDivergence;
 
    DataImport<number, dim> m_imCompressIndex;
 
    DataImport<MathVector<dim>, dim > m_imViscousForces[2];
 
 
  public:
      typedef SmartPtr<CplUserData<number, dim> > NumberExport;
      typedef SmartPtr<CplUserData<MathVector<dim>, dim> > VectorExport;
      typedef SmartPtr<CplUserData<MathMatrix<dim, dim>, dim> > MatrixExport;
 
      NumberExport divergence() {return m_exDivergence;}
      VectorExport displacement() {return m_exDisplacement;}
      MatrixExport stress() {return m_exStressTensor;}
 
 
  protected:
      SmartPtr<DataExport<number, dim> > m_exDivergence;
      SmartPtr<DataExport<MathVector<dim>, dim> >  m_exDisplacement;
      SmartPtr<DataExport<MathMatrix<dim,dim>, dim> >  m_exStressTensor;
 
 
      template <typename TElem, typename TFEGeom>
      void ex_divergence_fe(number vValue[],
          const MathVector<dim> vGlobIP[],
          number time, int si,
          const LocalVector& u,
          GridObject* elem,
          const MathVector<dim> vCornerCoords[],
          const MathVector<TFEGeom::dim> vLocIP[],
          const size_t nip,
          bool bDeriv,
          std::vector<std::vector<number> > vvvDeriv[]);
 
 
      template <typename TElem, typename TFEGeom>
      void ex_displacement_fe(MathVector<dim> vValue[],
          const MathVector<dim> vGlobIP[],
          number time, int si,
          const LocalVector& u,
          GridObject* elem,
          const MathVector<dim> vCornerCoords[],
          const MathVector<TFEGeom::dim> vLocIP[],
          const size_t nip,
          bool bDeriv,
          std::vector<std::vector<MathVector<dim> > > vvvDeriv[]);
 
      template <typename TElem, typename TFEGeom>
      void ex_stress_fe(MathMatrix<dim,dim> vValue[],
          const MathVector<dim> vGlobIP[],
          number time, int si,
          const LocalVector& u,
          GridObject* elem,
          const MathVector<dim> vCornerCoords[],
          const MathVector<TFEGeom::dim> vLocIP[],
          const size_t nip,
          bool bDeriv,
          std::vector<std::vector<MathMatrix<dim,dim> > > vvvDeriv[]);
 
 
};
 
// end group small_strain_mechanics
 
} //end of namespace SmallStrainMechanics
} //end of namespace ug
 
 
 
#endif /* SMALL_STRAIN_MECH_H_ */