nedelec_local_ass.h

Go to the documentation of this file.

/*
 * Copyright (c) 2013-2014:  G-CSC, Goethe University Frankfurt
 * Author: Dmitry Logashenko
 * 
 * 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.
 */
 
/*
 * nedelec_local_ass.h
 * Declarations of the functions for assembling the local stiffness, mass and
 * weak divergence matrices of the Nedelec-type-1 (Whitney-1) based discretization
 * of the Maxwell equations, as well as the shape functions and their curls
 * of the Nedelec element.
 */
 
#ifndef __H__UG__PLUGINS__ELECTROMAGNETISM__ROT_ROT_ASS__
#define __H__UG__PLUGINS__ELECTROMAGNETISM__ROT_ROT_ASS__
 
// other ug4 modules
#include "common/common.h"
#include "lib_grid/lg_base.h"
#include "lib_disc/reference_element/reference_element.h"
#include "lib_disc/reference_element/element_list_traits.h"
#include "lib_disc/local_finite_element/lagrange/lagrangep1.h"
 
namespace ug{
namespace Electromagnetism{
 
 
 
template <typename TDomain, typename TElem>
class NedelecT1_LDisc
{
public:
 
  static const int WDim = TDomain::dim;
  
  typedef typename TDomain::grid_type grid_type;
 
  typedef typename TDomain::position_type position_type;
  
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
  typedef LagrangeP1<ref_elem_type> W0_shapes_type;
 
  static const int dim = ref_elem_type::dim;
 
  static const size_t numCorners = ref_elem_type::numCorners;
  
  static const size_t numEdges = ref_elem_type::numEdges;
  
  static const size_t maxNumEdges = (size_t) element_list_traits<typename domain_traits<WDim>::DimElemList>::maxEdges;
 
public:
 
  static void local_stiffness_and_mass
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number S [maxNumEdges][maxNumEdges], 
    number M [maxNumEdges][maxNumEdges] 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
  
  static void local_mass
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number M [maxNumEdges][maxNumEdges] 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
  
  static void local_div_matrix
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number B [][numEdges] 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
  
  static void get_shapes
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const MathVector<dim> local, 
    MathVector<WDim> shapes [] 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
  
  static void interpolate
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const number dofs [], 
    const MathVector<dim> local [], 
    const size_t n_pnt, 
    MathVector<WDim> values [] 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
  
  static void curl
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const number dofs [], 
    MathVector<WDim> & curl_vec 
  )
  {
    UG_THROW ("Whitney-1 (Nedelec) shapes not implemented for roid " << ref_elem_type::REFERENCE_OBJECT_ID << ".");
  };
};
 
 
template <typename TDomain, typename TElem>
class NedelecT1_LDisc_forSimplex
{
public:
 
  static const int WDim = TDomain::dim;
  
  typedef typename TDomain::grid_type grid_type;
 
  typedef typename TDomain::position_type position_type;
  
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
  typedef LagrangeP1<ref_elem_type> W0_shapes_type;
 
  static const int dim = ref_elem_type::dim;
 
  static const size_t numCorners = ref_elem_type::numCorners;
  
  static const size_t numEdges = ref_elem_type::numEdges;
  
  static const size_t maxNumEdges = (size_t) element_list_traits<typename domain_traits<WDim>::DimElemList>::maxEdges;
 
private:
  
  static void compute_W0_grads
  (
    const position_type * corners, 
    MathVector<WDim> grad [numCorners] 
  );
  
  static void get_edge_corners
  (
    const TDomain * domain, 
    TElem * elem, 
    size_t edge_corner [numEdges] [2] 
  );
  
public:
  
  static void local_stiffness_and_mass
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number S [maxNumEdges][maxNumEdges], 
    number M [maxNumEdges][maxNumEdges] 
  );
  
  static void local_mass
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number M [maxNumEdges][maxNumEdges] 
  );
  
  static void local_div_matrix
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    number B [][numEdges] 
  );
  
  static void get_shapes
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const MathVector<dim> local, 
    MathVector<WDim> shapes [] 
  );
  
  static void interpolate
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const number dofs [], 
    const MathVector<dim> local [], 
    const size_t n_pnt, 
    MathVector<WDim> values [] 
  );
  
  static void curl
  (
    const TDomain * domain, 
    TElem * elem, 
    const position_type * corners, 
    const number dofs [], 
    MathVector<WDim> & curl_vec 
  );
};
 
 
template <typename TDomain>
class NedelecT1_LDisc<TDomain, Triangle> : public NedelecT1_LDisc_forSimplex<TDomain, Triangle>
{
};
 
 
template <typename TDomain>
class NedelecT1_LDisc<TDomain, Tetrahedron> : public NedelecT1_LDisc_forSimplex<TDomain, Tetrahedron>
{
};
 
 
 
template <typename TDomain, int refDim, int WDim = TDomain::dim>
class NedelecInterpolation
{
public:
 
  typedef typename TDomain::position_type position_type;
 
  static const size_t maxNumEdges = (size_t) element_list_traits<typename domain_traits<WDim>::DimElemList>::maxEdges;
 
public:
  static void value
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<refDim> local [], 
    const size_t n_pnt, 
    MathVector<WDim> values [] 
  )
  { // This is a generic version: It prints the error message. Cf. the specializations below.
    UG_THROW ("NedelecInterpolation::value: No implementation of the Nedelec elements for "
          "Reference Object " << elem->reference_object_id () <<
          " of reference dim. " << refDim << " in a " << WDim << "d domain.");
  }
  
 
  static void curl
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    MathVector<WDim> & curl_vec 
  )
  { // This is a generic version: It prints the error message. Cf. the specializations below.
    UG_THROW ("NedelecInterpolation::curl: No implementation of the Nedelec elements for "
          "Reference Object " << elem->reference_object_id () <<
          " of reference dim. " << refDim << " in a " << WDim << "d domain.");
  }
  
 
  static number curl_flux
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<WDim> & normal, 
    const position_type pnt, 
    number & flux 
  )
  { // This is a generic version: It prints the error message. Cf. the specializations below.
    UG_THROW ("NedelecInterpolation::curl_flux: No implementation of the Nedelec elements for "
          "Reference Object " << elem->reference_object_id () <<
          " of reference dim. " << refDim << " in a " << WDim << "d domain.");
    return false;
  }
};
 
template <typename TDomain>
class NedelecInterpolation<TDomain, 2, 2>
{
public:
 
  typedef typename TDomain::position_type position_type;
 
  static const size_t maxNumEdges = (size_t) element_list_traits<typename domain_traits<2>::DimElemList>::maxEdges;
 
public:
  static void value
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<2> local [], 
    const size_t n_pnt, 
    MathVector<2> values [] 
  )
  {
    if (elem->reference_object_id () != ROID_TRIANGLE)
      UG_THROW ("No implementation of the Nedelec elements for "
            "Reference Object " << elem->reference_object_id () <<
            " of reference dim. 2 in a 2d domain. (This must be a triangle.)");
    NedelecT1_LDisc<TDomain, Triangle>::interpolate
      (domain, (Triangle *) elem, corners, dofs, local, n_pnt, values);
  }
  
 
  static void curl
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    MathVector<2> & curl_vec 
  )
  {
    if (elem->reference_object_id () != ROID_TRIANGLE)
      UG_THROW ("No implementation of the Nedelec elements for "
            "Reference Object " << elem->reference_object_id () <<
            " of reference dim. 2 in a 2d domain. (This must be a triangle.)");
    NedelecT1_LDisc<TDomain, Triangle>::curl
      (domain, (Triangle *) elem, corners, dofs, curl_vec);
  }
  
 
  static number curl_flux
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<2> & normal, 
    const position_type pnt, 
    number & flux 
  )
  { // This is a generic version: It prints the error message. Cf. the specializations below.
    UG_THROW ("NedelecInterpolation::curl_flux: No implementation of the Nedelec elements for "
          "Reference Object " << elem->reference_object_id () <<
          " of reference dim. 3 in a 3d domain.");
    return false;
  }
};
 
template <typename TDomain>
class NedelecInterpolation<TDomain, 3, 3>
{
public:
 
  typedef typename TDomain::position_type position_type;
 
  static const size_t maxNumEdges = (size_t) element_list_traits<typename domain_traits<3>::DimElemList>::maxEdges;
 
public:
  static void value
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<3> local [], 
    const size_t n_pnt, 
    MathVector<3> values [] 
  )
  {
    if (elem->reference_object_id () != ROID_TETRAHEDRON)
      UG_THROW ("No implementation of the Nedelec elements for "
            "Reference Object " << elem->reference_object_id () <<
            " of reference dim. 3 in a 3d domain. (This must be a tetrahedron.)");
    NedelecT1_LDisc<TDomain, Tetrahedron>::interpolate
      (domain, (Tetrahedron *) elem, corners, dofs, local, n_pnt, values);
  }
  
 
  static void curl
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    MathVector<3> & curl_vec 
  )
  {
    if (elem->reference_object_id () != ROID_TETRAHEDRON)
      UG_THROW ("No implementation of the Nedelec elements for "
            "Reference Object " << elem->reference_object_id () <<
            " of reference dim. 3 in a 3d domain. (This must be a tetrahedron.)");
    NedelecT1_LDisc<TDomain, Tetrahedron>::curl
      (domain, (Tetrahedron *) elem, corners, dofs, curl_vec);
  }
  
 
  static number curl_flux
  (
    const TDomain * domain, 
    GridObject * elem, 
    const position_type corners [], 
    const number dofs [], 
    const MathVector<3> & normal, 
    const position_type pnt, 
    number & flux 
  );
};
 
template <int dim>
inline void get_ave_vector
(
  const size_t n_vec, 
  const MathVector<dim> * vec, 
  MathVector<dim> & ave 
)
{
  ave = 0;
  for (size_t i = 0; i < n_vec; i++) ave += vec [i];
  ave /= n_vec;
}
 
} // end namespace Electromagnetism
} // end namespace ug
 
#include "nedelec_local_ass_impl.h"
 
#endif /* __H__UG__PLUGINS__ELECTROMAGNETISM__ROT_ROT_ASS__ */
 
/* End of File */