save_tensor_data.cpp

Go to the documentation of this file.

/** \file save_tensor_data.cpp
 
  \brief Generate and save on tags thermal stress according to SSLV116 test
 
*/
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
inline auto voigt_notation() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'L', 9> L;
  FTensor::Dg<double, 3, 9> t_L;
  t_L(i, j, L) = 0;
  t_L(0, 0, 0) = 1;
  t_L(1, 1, 1) = 1;
  t_L(2, 2, 2) = 1;
  t_L(0, 1, 3) = 1;
  t_L(0, 2, 4) = 1;
  t_L(1, 2, 5) = 1;
  return t_L;
}
 
inline auto tensor_notation() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'L', 9> L;
  FTensor::Dg<double, 3, 9> t_L;
  t_L(i, j, L) = 0;
  t_L(0, 0, 0) = 1;
  t_L(0, 1, 1) = 1;
  t_L(0, 2, 2) = 1;
  t_L(1, 0, 3) = 1;
  t_L(1, 1, 4) = 1;
  t_L(1, 2, 5) = 1;
  t_L(2, 0, 6) = 1;
  t_L(2, 1, 7) = 1;
  t_L(2, 2, 8) = 1;
  return t_L;
}
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    // Create MoFEM database
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    char mesh_file[255] = "mesh.h5m";
    char output_file[255] = "out.h5m";
    char tag_name[255] = "INTERNAL_STRESS";
 
    int tag_len = 9;
 
    PetscBool use_voigt_notation = PETSC_FALSE;
 
    PetscOptionsBegin(m_field.get_comm(), "", "save thermal stress test",
                             "none");
    CHKERR PetscOptionsString("-file_name", " mesh file name", "", "mesh.h5m",
                              mesh_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-output_file", "output file name", "", "out.h5m",
                              output_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-tag_name", "Stress tag name", "",
                              "INTERNAL_STRESS", tag_name, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-voigt_notation", "use voigt notation", "",
                            use_voigt_notation, &use_voigt_notation,
                            PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG_TAG("WORLD", "save_tensor_data");
 
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    auto moab_comm_wrap =
        boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
 
    const std::string options = "PARALLEL=READ_PART;"
                                "PARALLEL_RESOLVE_SHARED_ENTS;"
                                "PARTITION=PARALLEL_PARTITION;";
    CHKERR moab.load_file(mesh_file, 0, options.c_str());
 
    Range src_elems;
    CHKERR pcomm->get_part_entities(src_elems, 3);
 
    Tag data_tag;
    double def_val[9];
    bzero(def_val, 9 * sizeof(double));
    CHKERR moab.tag_get_handle(tag_name, tag_len, MB_TYPE_DOUBLE, data_tag,
                               MB_TAG_CREAT | MB_TAG_DENSE, &def_val);
 
    // Set data on the mesh for SSLV116 test
    std::vector<double> spos;
    spos.resize(3 * src_elems.size());
    CHKERR moab.get_coords(src_elems, &spos[0]);
 
    VectorDouble data(9 * src_elems.size(), 0.0);
    auto t_data = getFTensor1FromArray<9, 9>(data);
 
    // temperature distribution for SSLV116 test
    // see https://code-aster.org/V2/doc/v14/en/man_v/v3/v3.04.116.pdf
    auto temperature = [](double z) -> double {
      return (z < -1)   ? (120.0 / -9.0) * (z + 10) + 250
             : (z <= 2) ? (50.0 / -3.0) * (z + 1) + 130
                        : (170.0 / 8.0) * (z - 2) + 80;
    };
 
    // material properties for SSLV116 test
    double young_modulus = 2e11;
    double poisson_ratio = 0.3;
    double alpha = 1e-5;
    double temp_0 = 250;
 
    double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
    double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
 
    // computing elasticity matrix
    MatrixDouble mat_D_ptr;
    mat_D_ptr.resize(6, 6, false);
    auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*(mat_D_ptr.data().begin()));
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
                      (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
                          t_kd(i, j) * t_kd(k, l);
 
    FTensor::Index<'L', 9> L;
    FTensor::Dg<double, 3, 9> t_L;
    if (use_voigt_notation) {
      t_L = voigt_notation();
    } else {
      t_L = tensor_notation();
    }
 
    FTensor::Tensor2<double, 3, 3> t_stress;
 
    bool tets_only = true;
 
    // computing thermal stress field (const on each element)
    int ielem = 0;
    for (auto sit = src_elems.begin(); sit != src_elems.end(); sit++) {
      double z = spos[3 * ielem + 2];
      double temp = temperature(z);
      t_stress(i, j) = -t_D(i, j, k, l) * t_kd(k, l) * (temp - temp_0) * alpha;
      t_data(L) = t_L(i, j, L) * t_stress(i, j);
      ++t_data;
      ++ielem;
      if (tets_only && type_from_handle(*sit) != MBTET) {
        tets_only = false;
      }
    }
 
    CHKERR moab.tag_set_data(data_tag, src_elems, &data[0]);
 
    // Check the integrated stress for the SSLV116 test
    if (tets_only) {
      std::vector<double> data_integ(tag_len, 0.0);
      for (auto &tet : src_elems) {
 
        std::vector<double> tet_data(tag_len, 0.0);
        CHKERR moab.tag_get_data(data_tag, &tet, 1, &tet_data[0]);
 
        const EntityHandle *vert_conn;
        int vert_num;
        CHKERR moab.get_connectivity(tet, vert_conn, vert_num, true);
        std::vector<double> vpos(3 * vert_num);
        CHKERR moab.get_coords(vert_conn, vert_num, &vpos[0]);
        double vol = Tools::tetVolume(vpos.data());
 
        for (int itag = 0; itag < tag_len; itag++) {
          data_integ[itag] += tet_data[itag] * vol;
        }
      }
 
      std::vector<double> global_data_integ(tag_len, 0.0);
      MPI_Allreduce(&data_integ[0], &global_data_integ[0], tag_len, MPI_DOUBLE,
                    MPI_SUM, m_field.get_comm());
 
      std::string temp = "Integrated stress for sslv116 test: ";
      for (int i = 0; i < 9; ++i) {
        std::stringstream ss;
        ss << std::scientific << std::setprecision(3) << global_data_integ[i];
        temp += ss.str() + " ";
      }
      MOFEM_LOG("WORLD", Sev::inform) << temp;
 
      double non_zero_val = 1.655e12;
      double non_zero_tol = 2e-3;
      std::vector<int> non_zero_inds(3);
      std::vector<int> zero_inds(6);
 
      if (use_voigt_notation) {
        non_zero_inds = {0, 1, 2};
        zero_inds = {3, 4, 5, 6, 7, 8};
      } else {
        non_zero_inds = {0, 4, 8};
         zero_inds = {1, 2, 3, 5, 6, 7};
      }
      bool non_zero_check =
          all_of(begin(non_zero_inds), end(non_zero_inds), [&](int i) {
            return abs(global_data_integ[i] - non_zero_val) / non_zero_val <
                   non_zero_tol;
          });
      bool zero_check = all_of(begin(zero_inds), end(zero_inds), [&](int i) {
        return abs(global_data_integ[i]) < 1e-12;
      });
      if (!non_zero_check || !zero_check) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Wrong value of the integrated stress");
      }
    } else {
      MOFEM_LOG("WORLD", Sev::warning)
          << "Mesh has non-tetrahedral elements, skipping the test based on "
             "integrating the stress in the volume";
    }
 
    CHKERR moab.write_file(output_file, "MOAB", "PARALLEL=WRITE_PART");
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}