assembly error with mixed meshes

[Andriy Andreykiv]

Good day!!

I have a model that is based on a mixed mesh of beams and plates. I am trying to assemble my Timoshenko beam torsional contribution on the beam region of the mesh.

When the mesh was uniform, with beams only, it worked fine, but when I added the plates (even though I assemble only the beam region) I'm starting to get dimension mismatches like this one:

Trace 1: Imported mesh with 128 beam and 96 plate elements
Def TangentProjectionGrad_Test_rotation:=(Grad_Test_rotation . Normalized(element_K)) . Normalized(element_K)
--------------------------------------------------------------------------------------^
Dot product of expressions of different sizes (2 != 3).

Kind request to advise.

Below I created a minimal working code with two GMSH meshes (beams only and a mixed mesh with beams and plates).

#include <catch2/catch.hpp>
 
#include <getfem/getfem_import.h>
#include <getfem/getfem_models.h>
 
std::string tensorTangentProjectionOnFrame(getfem::ga_workspace &workspace, const std::string &name)
{
  auto tangentName = "TangentProjection" + name;
  workspace.add_macro(tangentName, "(" + name + " . Normalized(element_K)) . Normalized(element_K)");
  return tangentName;
}
 
TEST_CASE("assembly with mixed meshes")
{
  auto meshName = "tower_with_floors.msh"; //mixed beam/plates mesh, doesn't work
  //auto meshName = "tower_dense.msh"; //beams only, works
  getfem::mesh mesh;
  std::map<std::string, bgeot::size_type> regionMap;
  std::set<bgeot::size_type> lowerDimConvexes;
  getfem::import_mesh_gmsh(meshName, mesh, true, &lowerDimConvexes,
                           &regionMap, false, nullptr, false);
 
  auto shellRegion = getfem::mesh_region::free_region_id(mesh);
  auto frameRegion = shellRegion + 1;
 
  size_t nFrames = 0;
  size_t nShells = 0;
 
  for (auto &&cv : dal::bv_iterable_c{mesh.convex_index()})
  {
    const auto pGeoTrans = mesh.trans_of_convex(cv);
    if (pGeoTrans->dim() == 1)
    {
      mesh.region(frameRegion).add(cv);
      ++nFrames;
    }
    else if (pGeoTrans->dim() == 2)
    {
      mesh.region(shellRegion).add(cv);
      ++nShells;
    }
  }
 
  GMM_SIMPLE_TRACE1("Imported mesh with " << nFrames << " beam and "
                                          << nShells << " plate elements");
  getfem::mesh_fem mf(mesh);
  mf.set_qdim(3);
  mf.set_classical_finite_element(2);
 
  getfem::mesh_im mim(mesh);
  mim.set_integration_method(4);
 
  getfem::model model;
  model.add_fem_variable("rotation", mf);
  model.add_initialized_scalar_data("G", 10e6); //shear modulus
  model.add_initialized_scalar_data("It", 1);   //torsional moment of inertia
 
  //Assembling Timoshenko torsional component on a frame region
  getfem::ga_workspace workspace(model);
  workspace.add_macro("Torsion", tensorTangentProjectionOnFrame(workspace, "Grad_rotation"));
  workspace.add_macro("VarTorsion", tensorTangentProjectionOnFrame(workspace, "Grad_Test_rotation"));
  workspace.add_expression("G * It * VarTorsion * Torsion", mim, frameRegion);
  
  REQUIRE_NOTHROW(workspace.assembly(1));
}

Best regards,

                        Andriy

tower_with_floors.msh
Description: Binary data

tower_dense.msh
Description: Binary data