map_disp_prism.cpp

Go to the documentation of this file.

/* Copyright (c) 2025 Authors under the MIT License.
 * See LICENSE.md for details.
 * SPDX-License-Identifier: MIT
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <iostream>
#include <fstream> // fstream
#include <vector>
#include <string>
#include <algorithm> // copy
#include <iterator>  // ostream_operator
#include <boost/tokenizer.hpp>
#include <DispMap.hpp>
 
using namespace boost;
using namespace CellEngineering;
 
static char help[] = "\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;
 
    // global variables
    char mesh_file_name[255];
    char data_file_name1[255];
    char data_file_name2[255];
    PetscBool flg_file = PETSC_FALSE;
    PetscInt order = 1;
    PetscBool flg_n_part = PETSC_FALSE;
    PetscInt n_partas = 1;
    PetscBool flg_thickness = PETSC_FALSE;
    PetscBool flg_legacy_thickness = PETSC_FALSE;
    double thickness = 0;
    double legacy_thickness = 0;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
    CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_file);
 
    CHKERR PetscOptionsString("-my_data_x", "data file name", "", "data_x.h5m",
                              data_file_name1, 255, PETSC_NULLPTR);
 
    CHKERR PetscOptionsString("-my_data_y", "data file name", "", "data_y.h5m",
                              data_file_name2, 255, PETSC_NULLPTR);
 
    CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 1,
                           &order, PETSC_NULLPTR);
 
    CHKERR PetscOptionsInt("-my_nparts", "number of parts", "", 1, &n_partas,
                           &flg_n_part);
 
    CHKERR PetscOptionsScalar("-my_thickness", "top layer thickness", "",
                              thickness, &thickness, &flg_thickness);
    CHKERR PetscOptionsScalar(
        "-my_thinckness", "deprecated alias for -my_thickness", "",
        legacy_thickness, &legacy_thickness, &flg_legacy_thickness);
 
    PetscOptionsEnd();
 
    if (flg_thickness != PETSC_TRUE &&
        flg_legacy_thickness == PETSC_TRUE) {
      thickness = legacy_thickness;
    }
 
    if (flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
    }
 
    if (flg_n_part != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR partitioning number not given");
    }
 
    if (flg_thickness != PETSC_TRUE &&
        flg_legacy_thickness != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR thickness of top layer not given");
    }
 
    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == nullptr)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
    // meshset consisting all entities in mesh
    EntityHandle root_set = moab.get_root_set();
 
    // Seed all mesh entities to MoFEM database, those entities can be
    // potentially used as finite elements or as entities which carry some
    // approximation field.
    BitRefLevel bit_level0;
    bit_level0.set(0);
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        root_set, 3, bit_level0);
    // Add field
    CHKERR m_field.add_field("DISP_X", H1, AINSWORTH_LEGENDRE_BASE, 1);
    CHKERR m_field.add_field("DISP_Y", H1, AINSWORTH_LEGENDRE_BASE, 1);
 
    auto *meshsets_manager_ptr = m_field.getInterface<MeshsetsManager>();
 
    Range surface;
    CHKERR meshsets_manager_ptr->getEntitiesByDimension(101, SIDESET, 2,
                                                        surface, true);
    CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_X");
    CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_Y");
 
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_X", 1);
    CHKERR m_field.set_field_order(0, MBEDGE, "DISP_X", order);
    CHKERR m_field.set_field_order(0, MBTRI, "DISP_X", order);
 
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_Y", 1);
    CHKERR m_field.set_field_order(0, MBEDGE, "DISP_Y", order);
    CHKERR m_field.set_field_order(0, MBTRI, "DISP_Y", order);
 
    CHKERR m_field.build_fields();
 
    CHKERR m_field.add_finite_element("DISP_MAP_ELEMENT");
    CHKERR m_field.modify_finite_element_add_field_row("DISP_MAP_ELEMENT",
                                                       "DISP_X");
    CHKERR m_field.modify_finite_element_add_field_col("DISP_MAP_ELEMENT",
                                                       "DISP_X");
    CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
                                                        "DISP_X");
    CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
                                                        "DISP_Y");
 
    // Add entities to that element
    CHKERR m_field.add_ents_to_finite_element_by_type(surface, MBTRI,
                                                      "DISP_MAP_ELEMENT");
 
    // build finite elements
    CHKERR m_field.build_finite_elements();
    // build adjacencies between elements and degrees of freedom
    CHKERR m_field.build_adjacencies(bit_level0);
 
    // set discrete manager
    DM dm_disp;
    DMType dm_name = "DM_DISP";
    // Register DM problem
    CHKERR DMRegister_MoFEM(dm_name);
    CHKERR DMCreate(PETSC_COMM_WORLD, &dm_disp);
    CHKERR DMSetType(dm_disp, dm_name);
    // Create DM instance
    CHKERR DMMoFEMCreateMoFEM(dm_disp, &m_field, dm_name, bit_level0);
    // Configure DM form line command options (DM itself, solvers,
    // pre-conditioners, ... )
    CHKERR DMSetFromOptions(dm_disp);
    // Add elements to dm (only one here)
    CHKERR DMMoFEMAddElement(dm_disp, "DISP_MAP_ELEMENT");
    // Set up problem (number dofs, partition mesh, etc.)
    CHKERR DMSetUp(dm_disp);
 
    DispMapFe fe_face(m_field);
 
    // string file(argv[1]); cout << "argument 1:" << argv[1] <<endl;
    // string file2(argv[2]); cout << "argument 2:" << argv[2] <<endl;
 
    string file1(data_file_name1);
    DataFromFiles data_from_files(file1);
    CHKERR data_from_files.fromOptions();
    CHKERR data_from_files.loadFileData();
 
    string file2(data_file_name2);
    DataFromFiles data_from_files2(file2);
    CHKERR data_from_files2.fromOptions();
    CHKERR data_from_files2.loadFileData();
 
    auto &list_of_operators = fe_face.getOpPtrVector();
 
    Mat A;
    Vec Fx, Dx;
    CHKERR DMCreateMatrix(dm_disp, &A);
    CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DM_DISP", ROW,
                                                              &Fx);
    CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DM_DISP", COL,
                                                              &Dx);
    Vec Fy, Dy;
    CHKERR VecDuplicate(Fx, &Fy);
    CHKERR VecDuplicate(Dx, &Dy);
 
    CHKERR MatZeroEntries(A);
    CHKERR VecZeroEntries(Fx);
    CHKERR VecZeroEntries(Fy);
 
    CommonData common_data;
    common_data.globA = A;
 
    common_data.globF = Fx;
 
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    list_of_operators.push_back(new OpCalculateHOJacForFace(jac_ptr));
    list_of_operators.push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    list_of_operators.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
    list_of_operators.push_back(
        new OpCalculateLhs("DISP_X", "DISP_X", common_data, data_from_files));
    list_of_operators.push_back(
        new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files));
    list_of_operators.push_back(new OpAssembleRhs("DISP_X", common_data));
 
    CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
 
    // cout << "error.no !" <<error++ << endl;
    CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    CHKERR VecGhostUpdateBegin(Fx, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(Fx, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(Fx);
    CHKERR VecAssemblyEnd(Fx);
 
    // int error=1;
    common_data.globF = Fy;
    list_of_operators.clear();
    list_of_operators.push_back(new OpCalculateHOJacForFace(jac_ptr));
    list_of_operators.push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    list_of_operators.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
    list_of_operators.push_back(
        new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files2));
    list_of_operators.push_back(new OpAssembleRhs("DISP_X", common_data));
    // cout << "error.no !"<< error++<< endl;
 
    CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
    // cout << "error.no !" <<error++ << endl;
    CHKERR VecGhostUpdateBegin(Fy, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(Fy, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(Fy);
    CHKERR VecAssemblyEnd(Fy);
 
    // VecView(F,PETSC_VIEWER_DRAW_WORLD);
    // MatView(A,PETSC_VIEWER_DRAW_WORLD);
    // std::string wait;
    // std::cin >> wait;
 
    // solve
    {
      KSP ksp;
      CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
      CHKERR KSPSetOperators(ksp, A, A);
      CHKERR KSPSetFromOptions(ksp);
      CHKERR KSPSolve(ksp, Fx, Dx);
      CHKERR KSPSolve(ksp, Fy, Dy);
    }
 
    CHKERR VecGhostUpdateBegin(Dx, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(Dx, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
        "DM_DISP", COL, Dx, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateBegin(Dy, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(Dy, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
        "DM_DISP", "DISP_X", "DISP_Y", COL, Dy, INSERT_VALUES, SCATTER_REVERSE);
 
    // Remove rigid translation
    double sum_dx = 0;
    int nb_dx = 0;
    for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
                                                       MBVERTEX, dof)) {
      sum_dx += (*dof)->getFieldData();
      nb_dx++;
    }
    sum_dx /= nb_dx;
    for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
                                                       MBVERTEX, dof)) {
      (*dof)->getFieldData() -= sum_dx;
    }
    double sum_dy = 0;
    int nb_dy = 0;
    for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
                                                       MBVERTEX, dof)) {
      sum_dy += (*dof)->getFieldData();
      nb_dy++;
    }
    sum_dy /= nb_dy;
    for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
                                                       MBVERTEX, dof)) {
      (*dof)->getFieldData() -= sum_dy;
    }
 
    using PostProcFace =
        PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;
    using OpPPMap = OpPostProcMapInMoab<2, 2>;
 
    PostProcFace post_proc(m_field);
    CHKERR AddHOOps<2, 2, 3>::add(post_proc.getOpPtrVector(), {H1});
 
    auto disp_x_ptr = boost::make_shared<VectorDouble>();
    auto disp_y_ptr = boost::make_shared<VectorDouble>();
    post_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("DISP_X", disp_x_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("DISP_Y", disp_y_ptr));
    post_proc.getOpPtrVector().push_back(new OpPPMap(
        post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
        OpPPMap::DataMapVec{{"DISP_X", disp_x_ptr},
                            {"DISP_Y", disp_y_ptr}},
        OpPPMap::DataMapMat{}, OpPPMap::DataMapMat{},
        OpPPMap::DataMapMat{}));
 
    CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &post_proc);
    CHKERR post_proc.writeFile("out_disp.h5m");
 
    CHKERR VecDestroy(&Fx);
    CHKERR VecDestroy(&Fy);
    CHKERR VecDestroy(&Dx);
    CHKERR VecDestroy(&Dy);
    CHKERR MatDestroy(&A);
    CHKERR DMDestroy(&dm_disp);
 
    // Delete element which will be no longer used
    CHKERR m_field.delete_finite_element("DISP_MAP_ELEMENT");
 
    // Add top layer
    {
      auto mmanager_ptr = m_field.getInterface<MeshsetsManager>();
      Range ents_2nd_layer;
      CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
                                                  ents_2nd_layer, true);
      Range tris_nodes;
      CHKERR moab.get_connectivity(ents_2nd_layer, tris_nodes, true);
      std::map<EntityHandle, EntityHandle> map_nodes;
      for (Range::iterator nit = tris_nodes.begin(); nit != tris_nodes.end();
           nit++) {
        double coords[3];
        CHKERR moab.get_coords(&*nit, 1, coords);
        coords[2] += thickness;
        CHKERR moab.create_vertex(coords, map_nodes[*nit]);
      }
 
      Range prisms;
      Range top_tris;
      for (Range::iterator tit = ents_2nd_layer.begin();
           tit != ents_2nd_layer.end(); tit++) {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR moab.get_connectivity(*tit, conn, num_nodes);
        EntityHandle conn_prism[6];
        for (int n = 0; n != 3; n++) {
          conn_prism[n] = conn[n];
          conn_prism[n + 3] = map_nodes[conn[n]];
        }
        EntityHandle prism;
        CHKERR moab.create_element(MBPRISM, conn_prism, 6, prism);
        Range adj;
        CHKERR moab.get_adjacencies(&prism, 1, 2, true, adj);
        CHKERR moab.get_adjacencies(&prism, 1, 1, true, adj);
        prisms.insert(prism);
        EntityHandle tri;
        CHKERR moab.side_element(prism, 2, 4, tri);
        top_tris.insert(tri);
      }
      CHKERR mmanager_ptr->addMeshset(SIDESET, 202);
      CHKERR mmanager_ptr->addEntitiesToMeshset(SIDESET, 202, top_tris);
      CHKERR mmanager_ptr->addMeshset(BLOCKSET, 2, "MAT_ELASTIC_2");
      CHKERR mmanager_ptr->addEntitiesToMeshset(BLOCKSET, 2, prisms);
      Mat_Elastic mat_elastic;
      mat_elastic.data.Young = 1;
      mat_elastic.data.Poisson = 0;
      CHKERR mmanager_ptr->setAttributesByDataStructure(BLOCKSET, 2,
                                                        mat_elastic);
    }
 
    {
      // Add one node to fix
      auto mmanager_ptr = m_field.getInterface<MeshsetsManager>();
      Range ents_1st_layer;
      CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
                                                  ents_1st_layer, true);
      int num_nodes;
      const EntityHandle *conn;
      CHKERR moab.get_connectivity(ents_1st_layer[0], conn, num_nodes);
      Skinner skin(&moab);
      Range skin_edges;
      CHKERR skin.find_skin(0, ents_1st_layer, false, skin_edges);
      Range edges;
      CHKERR moab.get_adjacencies(&*ents_1st_layer.begin(), 1, 1, false, edges);
      EntityHandle meshset;
      CHKERR mmanager_ptr->getMeshset(202, SIDESET, meshset);
      CHKERR moab.add_entities(meshset, conn, 1);
      CHKERR moab.add_entities(meshset, skin_edges);
    }
 
    {
      Range ents3d;
      CHKERR moab.get_entities_by_dimension(0, 3, ents3d, false);
      CHKERR m_field.getInterface<CommInterface>()->partitionMesh(
          ents3d, 3, 2, n_partas);
    }
 
    CHKERR moab.write_file("analysis_mesh.h5m");
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  PetscFunctionReturn(0);
}