create_mwls_stress_mesh.cpp

Go to the documentation of this file.

/* 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/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
static char help[] = "Create mesh for MWLS approximation for internal stresses and density"
                     "\n\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, PETSC_NULL, help);
 
  try {
 
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
                                 255, &flg);
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);
 
    std::string mwls_stress_tag_name = "MED_SIGP";
    std::string mwls_rho_tag_name = "RHO_CELL";
    Tag th, th_rho;
    double def[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0};
    CHKERR moab.tag_get_handle(mwls_stress_tag_name.c_str(), 9, MB_TYPE_DOUBLE,
                               th, MB_TAG_CREAT | MB_TAG_SPARSE, def);
    CHKERR moab.tag_get_handle(mwls_rho_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
                               th_rho, MB_TAG_CREAT | MB_TAG_SPARSE, def);
 
    Range three_dimensional_entities;
    CHKERR moab.get_entities_by_dimension(0, 3, three_dimensional_entities,
                                          true);
    MatrixDouble coords;
    for (auto ent : three_dimensional_entities) {
      const EntityHandle *conn;
      int num_nodes;
      CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
      coords.resize(num_nodes, 3, false);
      CHKERR moab.get_coords(conn, num_nodes, &*coords.data().begin());
      auto get_centroid = [](auto &coords) {
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords(
            &coords(0, 0), &coords(0, 1), &coords(0, 2));
        FTensor::Index<'i', 3> i;
        FTensor::Tensor1<double, 3> t_center;
        t_center(i) = 0;
        for (int rr = 0; rr != coords.size1(); ++rr) {
          t_center(i) += t_coords(i);
          ++t_coords;
        }
        t_center(i) /= coords.size1();
        return t_center;
      };
      auto t_centre = get_centroid(coords);
 
      auto get_stress = [](auto t_center) {
        VectorDouble voight_stress(9);
        voight_stress(0) = 0.5*t_center(1); //t_center(0) * t_center(0) * t_center(0);
        voight_stress(1) = 0.5*t_center(1); //t_center(1) * t_center(1) * t_center(1);
        voight_stress(2) = 0.5*t_center(1); //t_center(2) * t_center(2) + t_center(2);
        voight_stress(3) = 0; //(t_center(0) * t_center(1)) * t_center(2);
        voight_stress(4) = 0; //(t_center(0) * t_center(2)) * t_center(1);
        voight_stress(5) = 0; //(t_center(1) * t_center(2)) * t_center(0);
        voight_stress(6) = 0; // dummy only used to make tag visable on ParaView
        voight_stress(7) = 0; // dummy only used to make tag visable on ParaView
        voight_stress(8) = 0; // dummy only used to make tag visable on ParaView
        return voight_stress;
      };
 
      auto get_density = [](auto t_center) {
        double density = 0.;
        density = 1 + t_center(0) * 2;
        return density;
      };
 
      auto voight_stress = get_stress(t_centre);
      double density = get_density(t_centre);
 
      CHKERR moab.tag_set_data(th, &ent, 1, &*voight_stress.data().begin());
      CHKERR moab.tag_set_data(th_rho, &ent, 1, &density);
    }
 
    EntityHandle meshset = 0;
    CHKERR moab.write_file("out_for_mwls.h5m", "MOAB", "", &meshset, 1);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}