Export im_nodes scalars on boundary regions

Questions and issues
1

Hello.
If I have a vector of scalars inside an im_data named “im_nodes_scalars” defined at a boundary region,

md.add_im_data(“im_nodes_scalars”, bounday_region_number)

representing the boundary condition of my problem (for example it could be an heat flow calculated in a sophisticated way), how can I visualize it?
I tried to export the scalars to vtu, making before a local projection of the im_nodes values on a discoutinuous bounday mesh fem (mfout), like the following

md.local_projection(mim, “im_nodes_scalars”, mfout, boundary_region_number)

but without success. An error is risen as “Singular system”.
I’m not sure whether local_projection works with boundary im_nodes or not.
Is there any other way to proceed?
Thank you.

2

Hello.
At the end I projected the boundary im_nodes on the boundary fem nodes by mean of the Galerkin method, declaring a filtered fem variable on the boundary region and setting a weak form of my boundary integral, eventually solving the linear system.
Anyhow, I’ve not understood why local_projection didn’t work.
Thank you,
Lorenzo

3

Hello Kostas,

I would like to share a short example showing my issue with the followin error.

RuntimeError: (Getfem::InterfaceError) – Error in ./gmm/gmm_dense_lu.h, line 183 :
Singular system, pivot = 4

Do you know why it happens?
Is there any other way to export im_nodes, without having to generate a partial fem and solving it.

import getfem as gf

# generate mesh
mesh = gf.Mesh("import", "structured",
               f"GT='GT_QK(2,1)';ORG=[0,0];SIZES=[{1},{1}];NSUBDIV=[{3},{3}]")
BRG_id = 1
mesh.set_region(BRG_id, mesh.outer_faces())

# set fem
mfout = gf.MeshFem(mesh, 1) 
mfout.set_classical_discontinuous_fem(1)

# set mim
im_degree = 3
mim = gf.MeshIm(mesh, im_degree)
mimd_BRG = gf.MeshImData(mim, BRG_id)

# generate model and data
md = gf.Model("real")
md.add_im_data("BRG_im_scalars", mimd_BRG)

scalars = md.interpolation("Norm(X)", mimd_BRG, BRG_id)
md.set_variable("BRG_im_scalars", scalars)
print(md.variable("BRG_im_scalars"))

# project scalars on a discontinuos mesh for output
md.local_projection(mim, "BRG_im_scalars", mfout, BRG_id)
4

Hi Lorenzo

Good that you have found a workaround, your solution is good, I would do the same.

Regarding local projection, this is its implementation in the code:

  void ga_local_projection(const getfem::model &md, const mesh_im &mim,
                           const std::string &expr, const mesh_fem &mf,
                           base_vector &result, const mesh_region &region) {

    // The case where the expression is a vector one and mf a scalar fem is
    // not taken into account for the moment.

    // Could be improved by not performing the assembly of the global mass matrix
    // working locally. This means a specific assembly.
    size_type nbdof = mf.nb_dof();
    model_real_sparse_matrix M(nbdof, nbdof);
    asm_mass_matrix(M, mim, mf, region);
    // FIXME: M should be cached for performance

    ga_workspace workspace(md, ga_workspace::inherit::NONE);
    gmm::sub_interval I(0,nbdof);
    workspace.add_fem_variable("c__dummy_var_95_", mf, I, base_vector(nbdof));
    if (mf.get_qdims().size() > 1)
      workspace.add_expression("("+expr+"):Test_c__dummy_var_95_",mim,region,2);
    else
      workspace.add_expression("("+expr+").Test_c__dummy_var_95_",mim,region,2);
    getfem::base_vector F(nbdof);
    workspace.set_assembled_vector(F);
    workspace.assembly(1);
    gmm::resize(result, nbdof);

    getfem::base_matrix loc_M;
    getfem::base_vector loc_U;
    for (mr_visitor v(region, mf.linked_mesh(), true); !v.finished(); ++v) {
      if (mf.convex_index().is_in(v.cv())) {
        size_type nd = mf.nb_basic_dof_of_element(v.cv());
        loc_M.base_resize(nd, nd); gmm::resize(loc_U, nd);
        gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));
        gmm::copy(gmm::sub_matrix(M, J, J), loc_M);
        gmm::lu_solve(loc_M, loc_U, gmm::sub_vector(F, J));
        gmm::copy(loc_U, gmm::sub_vector(result, J));
      }
    }
    MPI_SUM_VECTOR(result);
  }

In the line

    gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));

we can see that only entire elements are treated and not element faces. We should extend the implementation to cover this use case.

5

Hello Kostas,
thank you for the clear explanation.
Maybe this use case could become a new feature in the future.

Just for information, I tried also to introduce a partial mesh on the boundary and make a local_projeciton on it, adding the following rows of code:

mfout_BRG = gf.MeshFem('partial', mfout, mfout.basic_dof_on_region(1))
md.local_projection(mim, "BRG_im_scalars", mfout_BRG)

but also in this case there’s a Runtime Error as follows:

Sorry, cannot apply to reduced fems

Thank you,
Lorenzo