DataStructures.cpp

Go to the documentation of this file.

/** \file DataStructures.cpp
\brief Implementation for Data Structures in Forces and Sources

*/

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */

#ifdef __cplusplus
extern "C" {
#endif
  #include <cblas.h>
  #include <lapack_wrap.h>
  #include <gm_rule.h>
#ifdef __cplusplus
}
#endif

namespace MoFEM {

DataForcesAndSourcesCore::EntData::EntData()
    : sEnse(0), oRder(0), bAse(NOBASE) {
  for (int b = 0; b != LASTBASE;++b) {
    N[b].reset(new MatrixDouble());
    diffN[b].reset(new MatrixDouble());
  }
}

static void constructor_data(DataForcesAndSourcesCore *data,
                             const EntityType type) {

  using EntData = DataForcesAndSourcesCore::EntData;

  data->dataOnEntities[MBENTITYSET].push_back(new EntData());

  switch (type) {
  case MBENTITYSET:
    break;
  case MBTET: 
    data->dataOnEntities[MBVERTEX].push_back(new EntData());
    for (int ee = 0; ee != 6; ++ee) {
      data->dataOnEntities[MBEDGE].push_back(new EntData());
    }
    for (int ff = 0; ff != 4; ++ff) {
      data->dataOnEntities[MBTRI].push_back(new EntData());
    }
    data->dataOnEntities[MBTET].push_back(new EntData());
  break;
  case MBTRI:
    data->dataOnEntities[MBVERTEX].push_back(new EntData());
    for (int ee = 0; ee != 3; ++ee) {
      data->dataOnEntities[MBEDGE].push_back(new EntData());
    }
    data->dataOnEntities[MBTRI].push_back(new EntData());
  break;
  case MBEDGE:
    data->dataOnEntities[MBVERTEX].push_back(new EntData());
    data->dataOnEntities[MBEDGE].push_back(new EntData());
    break;
  case MBVERTEX:
    data->dataOnEntities[MBVERTEX].push_back(new EntData());
    break;
  case MBPRISM:
    data->dataOnEntities[MBVERTEX].push_back(new EntData());
    for (int ee = 0; ee != 9; ++ee) {
      data->dataOnEntities[MBEDGE].push_back(new EntData());
    }
    for (int qq = 0; qq != 5; ++qq) {
      data->dataOnEntities[MBQUAD].push_back(new EntData());
    }
    for (int ff = 0; ff != 5; ++ff) {
      data->dataOnEntities[MBTRI].push_back(new EntData());
    }
    data->dataOnEntities[MBPRISM].push_back(new EntData());
  break;
  default:
    throw MoFEMException(MOFEM_NOT_IMPLEMENTED);
  }
}

DataForcesAndSourcesCore::DataForcesAndSourcesCore(EntityType type) {
  constructor_data(this, type);
}

MoFEMErrorCode DataForcesAndSourcesCore::setElementType(const EntityType type) {
  MoFEMFunctionBegin;
  for (EntityType tt = MBVERTEX; tt != MBMAXTYPE; ++tt)
    dataOnEntities[tt].clear();
  constructor_data(this, type);
  MoFEMFunctionReturn(0);
}

static void constructor_derived_data(
    DerivedDataForcesAndSourcesCore *derived_data,
    const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr) {

  using EntData = DataForcesAndSourcesCore::EntData;
  using DerivedEntData = DerivedDataForcesAndSourcesCore::DerivedEntData;

  for (int tt = MBVERTEX; tt != MBMAXTYPE; ++tt) {
    auto &ent_data = data_ptr->dataOnEntities[tt];
    auto &derived_ent_data = derived_data->dataOnEntities[tt];
    for (auto &e : ent_data) {
      boost::shared_ptr<EntData> ent_data_ptr(data_ptr, &e);
      derived_ent_data.push_back(new DerivedEntData(ent_data_ptr));
    }
  }
}

DerivedDataForcesAndSourcesCore::DerivedDataForcesAndSourcesCore(
    const boost::shared_ptr<DataForcesAndSourcesCore> &data_ptr)
    : DataForcesAndSourcesCore(), dataPtr(data_ptr) {
  constructor_derived_data(this, dataPtr);
}

MoFEMErrorCode
DerivedDataForcesAndSourcesCore::setElementType(const EntityType type) {
  MoFEMFunctionBegin;
  for (EntityType tt = MBVERTEX; tt != MBMAXTYPE; ++tt)
    dataOnEntities[tt].clear();
  constructor_derived_data(this, dataPtr);
  MoFEMFunctionReturn(0);
}

std::ostream &operator<<(std::ostream &os,
                         const DataForcesAndSourcesCore::EntData &e) {
  os << "sEnse: " << e.getSense() << std::endl
     << "oRder: " << e.getOrder() << std::endl
     << "global indices: " << e.getIndices() << std::endl
     << "local indices: " << e.getLocalIndices() << std::endl;
  os.precision(2);
  os << "fieldData: " << std::fixed << e.getFieldData() << std::endl;
  MatrixDouble base = e.getN();
  MatrixDouble diff_base = e.getDiffN();
  const double eps = 1e-6;
  for (unsigned int ii = 0; ii != base.size1(); ii++) {
    for (unsigned int jj = 0; jj != base.size2(); jj++) {
      if (fabs(base(ii, jj)) < eps)
        base(ii, jj) = 0;
    }
  }
  for (unsigned int ii = 0; ii != diff_base.size1(); ii++) {
    for (unsigned int jj = 0; jj != diff_base.size2(); jj++) {
      if (fabs(diff_base(ii, jj)) < eps)
        diff_base(ii, jj) = 0;
    }
  }
  os << "N: " << std::fixed << base << std::endl
     << "diffN: " << std::fixed << diff_base;
  return os;
}

std::ostream &operator<<(std::ostream &os, const DataForcesAndSourcesCore &e) {
  for (unsigned int nn = 0; nn < e.dataOnEntities[MBVERTEX].size(); nn++) {
    os << "dataOnEntities[MBVERTEX][" << nn << "]" << std::endl
       << e.dataOnEntities[MBVERTEX][nn] << std::endl;
  }
  for (unsigned int ee = 0; ee < e.dataOnEntities[MBEDGE].size(); ee++) {
    os << "dataOnEntities[MBEDGE][" << ee << "]" << std::endl
       << e.dataOnEntities[MBEDGE][ee] << std::endl;
  }
  for (unsigned int ff = 0; ff < e.dataOnEntities[MBTRI].size(); ff++) {
    os << "dataOnEntities[MBTRI][" << ff << "] " << std::endl
       << e.dataOnEntities[MBTRI][ff] << std::endl;
  }
  for (unsigned int vv = 0; vv < e.dataOnEntities[MBTET].size(); vv++) {
    os << "dataOnEntities[MBTET][" << vv << "]" << std::endl
       << e.dataOnEntities[MBTET][vv] << std::endl;
  }
  return os;
}

/** \name Specializations for H1/L2 */

/**@{*/

template <>
FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
DataForcesAndSourcesCore::EntData::getFTensor1FieldData<3>() {
  if (dOfs[0]->getNbOfCoeffs() != 3) {
    std::stringstream s;
    s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
    s << " but you ask for tensor rank 1 dimension 3";
    THROW_MESSAGE(s.str());
  }
  double *ptr = &*fieldData.data().begin();
  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
                                                            &ptr[2]);
}

template <>
FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>
DataForcesAndSourcesCore::EntData::getFTensor1FieldData<2>() {
  if (dOfs[0]->getNbOfCoeffs() != 2) {
    std::stringstream s;
    s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
    s << " but you ask for tensor rank 1 dimension 3";
    THROW_MESSAGE(s.str());
  }
  double *ptr = &*fieldData.data().begin();
  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(ptr, &ptr[1]);
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>
DataForcesAndSourcesCore::EntData::getFTensor2FieldData<3, 3>() {
  if (dOfs[0]->getNbOfCoeffs() != 9) {
    std::stringstream s;
    s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
    s << " but you ask for tensor rank 2 dimensions 3 by 3 so 9 coefficients "
         "is expected";
    THROW_MESSAGE(s.str());
  }
  double *ptr = &*fieldData.data().begin();
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
      ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
      &ptr[8]);
}

template <>
FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
DataForcesAndSourcesCore::EntData::getFTensor2SymmetricFieldData<3>() {
  if (dOfs[0]->getNbOfCoeffs() != 6) {
    std::stringstream s;
    s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
    s << " but you ask for symmetric tensor rank 2 dimensions 3 by 3 so 6 "
         "coefficients "
         "is expected";
    THROW_MESSAGE(s.str());
  }
  double *ptr = &*fieldData.data().begin();
  return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>(
      ptr, &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
}

FTensor::Tensor0<FTensor::PackPtr<double *,1> >
DataForcesAndSourcesCore::EntData::getFTensor0FieldData() {
  if (dOfs[0]->getNbOfCoeffs() != 1) {
    std::stringstream s;
    s << "Wrong number of coefficients is " << dOfs[0]->getNbOfCoeffs();
    s << " but expected scalar field, tensor of rank 0";
    THROW_MESSAGE(s.str());
  }
  return FTensor::Tensor0<FTensor::PackPtr<double *, 1> >(
      &*fieldData.data().begin());
}

template <int Tensor_Dim>
FTensor::Tensor1<double *, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN(
    const FieldApproximationBase base) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor1<double *, Tensor_Dim>();
}

template <int Tensor_Dim>
FTensor::Tensor1<double *, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN() {
  return getFTensor1DiffN<Tensor_Dim>(bAse);
}

template <int Tensor_Dim>
FTensor::Tensor1<double *, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN(
    const FieldApproximationBase base, const int bb) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor1<double *, Tensor_Dim>();
}

template <int Tensor_Dim>
FTensor::Tensor1<double *, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN(const int bb) {
  return getFTensor1DiffN<Tensor_Dim>(bAse, bb);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
    const FieldApproximationBase base) {
  double *ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>() {
  return getFTensor1DiffN<3>(bAse);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
    const FieldApproximationBase base, const int bb) {
  double *ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor1<double *, 3>(
      &ptr[3 * bb], &ptr[3 * bb + 1], &ptr[3 * bb + 2], getDiffN(base).size2());
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int bb) {
  return getFTensor1DiffN<3>(bAse, bb);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 2d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
    const FieldApproximationBase base) {
  double *ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 2d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>() {
  return getFTensor1DiffN<2>(bAse);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
    const FieldApproximationBase base, const int bb) {
  double *ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor1<double *, 2>(&ptr[2 * bb], &ptr[2 * bb + 1],
                                       getDiffN(base).size1());
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int bb) {
  return getFTensor1DiffN<2>(bAse, bb);
}

template <int Tensor_Dim>
FTensor::Tensor1<double *, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN(
    const FieldApproximationBase base, const int gg, const int bb) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor1<double *, Tensor_Dim>();
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(
    const FieldApproximationBase base, const int gg, const int bb) {
  double *ptr = &getDiffN(base)(gg, 3 * bb);
  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2], 3);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 3d
 */
template <>
FTensor::Tensor1<double *, 3>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<3>(const int gg,
                                                       const int bb) {
  return getFTensor1DiffN<3>(bAse, gg, bb);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 2d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(
    const FieldApproximationBase base, const int gg, const int bb) {
  double *ptr = &getDiffN(base)(gg, 2 * bb);
  return FTensor::Tensor1<double *, 2>(ptr, &ptr[1], 2);
}

/**
 * \brief Get spatial derivative of base function tensor for dimension 2d
 */
template <>
FTensor::Tensor1<double *, 2>
DataForcesAndSourcesCore::EntData::getFTensor1DiffN<2>(const int gg,
                                                       const int bb) {
  return getFTensor1DiffN<2>(bAse, gg, bb);
}

/**@}*/

/** \name Specializations for HDiv/HCrul */

/**@{*/

template <int Tensor_Dim>
FTensor::Tensor1<FTensor::PackPtr<double *, Tensor_Dim>, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1N(
    FieldApproximationBase base) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor1<FTensor::PackPtr<double *, Tensor_Dim>, Tensor_Dim>();
}

template <int Tensor_Dim>
FTensor::Tensor1<FTensor::PackPtr<double *, Tensor_Dim>, Tensor_Dim>
DataForcesAndSourcesCore::EntData::getFTensor1N(FieldApproximationBase base,
                                                    const int gg,
                                                    const int bb) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor1<FTensor::PackPtr<double *, Tensor_Dim>, Tensor_Dim>();
}

template <int Tensor_Dim0, int Tensor_Dim1>
FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                 Tensor_Dim0, Tensor_Dim1>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN(
    FieldApproximationBase base) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
    << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor2<double *, Tensor_Dim0, Tensor_Dim1>();
}

template <int Tensor_Dim0, int Tensor_Dim1>
FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                 Tensor_Dim0, Tensor_Dim1>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN(
    FieldApproximationBase base, const int gg, const int bb) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim0 << "x" << Tensor_Dim1
    << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor2<double *, Tensor_Dim0, Tensor_Dim1>();
}

template <>
FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
DataForcesAndSourcesCore::EntData::getFTensor1N<3>(
    FieldApproximationBase base) {
  double *t_n_ptr = &*getN(base).data().begin();
  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
                                                            &t_n_ptr[HVEC1],
                                                            &t_n_ptr[HVEC2]);
}

template <>
FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
DataForcesAndSourcesCore::EntData::getFTensor1N<3>(
    FieldApproximationBase base, const int gg, const int bb) {
  double *t_n_ptr = &getN(base)(gg, 3 * bb);
  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(t_n_ptr, // HVEC0
                                                            &t_n_ptr[HVEC1],
                                                            &t_n_ptr[HVEC2]);
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
    FieldApproximationBase base) {
  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
      t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
      &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
      &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 3>(
    FieldApproximationBase base, const int gg, const int bb) {
  double *t_diff_n_ptr = &getDiffN(base)(gg, 9 * bb);
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
      t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC0_2],
      &t_diff_n_ptr[HVEC1_0], &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC1_2],
      &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1], &t_diff_n_ptr[HVEC2_2]);
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
    FieldApproximationBase base) {
  double *t_diff_n_ptr = &*getDiffN(base).data().begin();
  return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>(
      t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
      &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0],
      &t_diff_n_ptr[HVEC2_1]);
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>
DataForcesAndSourcesCore::EntData::getFTensor2DiffN<3, 2>(
    FieldApproximationBase base, const int gg, const int bb) {
  double *t_diff_n_ptr = &getDiffN(base)(gg, 6 * bb);
  return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>(
      t_diff_n_ptr, &t_diff_n_ptr[HVEC0_1], &t_diff_n_ptr[HVEC1_0],
      &t_diff_n_ptr[HVEC1_1], &t_diff_n_ptr[HVEC2_0], &t_diff_n_ptr[HVEC2_1]);
}

template <int Tensor_Dim0, int Tensor_Dim1>
FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                 Tensor_Dim0, Tensor_Dim1>
DataForcesAndSourcesCore::EntData::getFTensor2N(FieldApproximationBase base) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
    << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                          Tensor_Dim0, Tensor_Dim1>();
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>
DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
    FieldApproximationBase base) {
  double *t_n_ptr = &*(getN(base).data().begin());
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(

      &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],

      &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],

      &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]

  );
}

template <int Tensor_Dim0, int Tensor_Dim1>
FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                 Tensor_Dim0, Tensor_Dim1>
DataForcesAndSourcesCore::EntData::getFTensor2N(FieldApproximationBase base,
                                                const int gg, const int bb) {
  std::stringstream s;
  s << "Template for tensor dimension " << Tensor_Dim0 << ", " << Tensor_Dim1
    << " not implemented";
  THROW_MESSAGE(s.str());
  return FTensor::Tensor2<FTensor::PackPtr<double *, Tensor_Dim0 * Tensor_Dim1>,
                          Tensor_Dim0, Tensor_Dim1>();
}

template <>
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>
DataForcesAndSourcesCore::EntData::getFTensor2N<3, 3>(
    FieldApproximationBase base, const int gg, const int bb) {
  double *t_n_ptr = &getN(base)(gg, 9 * bb);
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(

      &t_n_ptr[HVEC0], &t_n_ptr[HVEC1], &t_n_ptr[HVEC2],

      &t_n_ptr[3 + HVEC0], &t_n_ptr[3 + HVEC1], &t_n_ptr[3 + HVEC2],

      &t_n_ptr[6 + HVEC0], &t_n_ptr[6 + HVEC1], &t_n_ptr[6 + HVEC2]

  );
}

/**@}*/

} // namespace MoFEM