[Getfem-users] [GetFem-users] Problem in Dof Enumeration and consequent wrong matrix assembly

[Simone Di Gregorio]

Dear GetFem++ users,

I found a problem regarding the resolution of a fem problem. I must solve a 1D-network in terms of pressure and velocity. I have a 1-D network whose mesh structure, made with regions, is like below

BEGIN MESH STRUCTURE DESCRIPTION

CONVEX 0    'GT_PK(1,1)'      0  1
CONVEX 1    'GT_PK(1,1)'      1  2
CONVEX 2    'GT_PK(1,1)'      2  3
CONVEX 3    'GT_PK(1,1)'      4  5
CONVEX 4    'GT_PK(1,1)'      5  6
CONVEX 5    'GT_PK(1,1)'      6  3
CONVEX 6    'GT_PK(1,1)'      3  7
CONVEX 7    'GT_PK(1,1)'      7  8
CONVEX 8    'GT_PK(1,1)'      8  9

END MESH STRUCTURE DESCRIPTION
BEGIN REGION 0
0 1 2
END REGION 0
BEGIN REGION 1
3 4 5
END REGION 1
BEGIN REGION 2
6 7 8
END REGION 2

I use this part of code to obtain mesh_fem:

// pf_Uvi is 'FEM_PK(1,1)' ;

//pf_coefvi is 'FEM_PK(1,0)';

//meshv is the mesh described before;
//nb_branches=3;

        for(size_type i=0; i<nb_branches; ++i){

        mesh_fem mf_tmp(meshv);
        mf_tmp.set_finite_element(meshv.region(i).index(), pf_coefv);
        mf_coefvi.emplace_back(mf_tmp);
        mf_tmp.clear();

        mesh_fem mf_tmp1(meshv);
        mf_tmp1.set_finite_element(meshv.region(i).index(), pf_Uv);
        mf_Uvi.emplace_back(mf_tmp1);
        mf_tmp1.clear();
    }

It usually works well for various meshes but in this case in the REGION 1 the dof enumeration of mf_Uvi is wrong, in fact  it comes out like:

BEGIN MESH_FEM

QDIM 1
 CONVEX 3 'FEM_PK(1,1)'
 CONVEX 4 'FEM_PK(1,1)'
 CONVEX 5 'FEM_PK(1,1)'
 BEGIN DOF_ENUMERATION
  3:  3 2
  4:  2 0
  5:  0 1
 END DOF_ENUMERATION
END MESH_FEM

So the Dof enumeration is 3 2   2 0   0 1 while in all other branches is 0 1   1 2   2 3.

While for mf_coefvi in REGION 1 the dof enumeration is in descending order instead of being in ascending order. (So it is 2 1 0 instead of 0 1 2)

BEGIN MESH_FEM

QDIM 1
 CONVEX 3 'FEM_PK(1,0)'
 CONVEX 4 'FEM_PK(1,0)'
 CONVEX 5 'FEM_PK(1,0)'
 BEGIN DOF_ENUMERATION
  3:  2
  4:  1
  5:  0
 END DOF_ENUMERATION
END MESH_FEM

The problem is not only here, but also if I use this mesh_fem to build the mass matrix, it results to be not symmetric.
 
//descr.IM_TYPEV='IM_GAUSS1D(6)'

    pintegration_method pim_v = int_method_descriptor(descr.IM_TYPEV);
    mimv.set_integration_method(meshv.convex_index(), pim_v);

    getfem::asm_mass_matrix_param(M, mimv, mf_Uvi[1], mf_coefvi[1], coef, meshv.region(1));

Moreover, also the results of the fem problem are completely wrong.
Instead I noticed that modifying the convex of REGION 1 of the mesh in this way (inverting the node numbers):

CONVEX 3    'GT_PK(1,1)'      3  6
CONVEX 4    'GT_PK(1,1)'      6  5
CONVEX 5    'GT_PK(1,1)'      5  4

I obtain:
  - a correct Dof Enumeration:
BEGIN DOF_ENUMERATION
  3:  0 1
  4:  1 2
  5:  2 3
 END DOF_ENUMERATION
- a symmetric mass matrix
- a correct solution in terms of velocity
but a wrong solution in terms of pressure (maybe for the boundary conditions).
So I would to know if it is possible to solve the problem and what could be the reason of this error.

Best regards,

Simone Di Gregorio