Problem during boundary condition assembly

Questions and issues
101

Hello Kostas,

Got it - I am looking into that.

By the way, my regions are defined as

# detect the border of the mesh
fb0             = m.outer_faces_with_direction([ 1., 0., 0., 0., 0., 0.], 0.01)
fb1             = m.outer_faces_with_direction([-1., 0., 0., 0., 0., 0.], 0.01)
fb2             = m.outer_faces_with_direction([0.,  1., 0., 0., 0., 0.], 0.01)
fb3             = m.outer_faces_with_direction([0., -1., 0., 0., 0., 0.], 0.01)
fb4             = m.outer_faces_with_direction([0., 0.,  1., 0., 0., 0.], 0.01)
fb5             = m.outer_faces_with_direction([0., 0., -1., 0., 0., 0.], 0.01)
fb6             = m.outer_faces_with_direction([0., 0., 0.,  1., 0., 0.], 0.01)
fb7             = m.outer_faces_with_direction([0., 0., 0., -1., 0., 0.], 0.01)
fb8             = m.outer_faces_with_direction([0., 0., 0., 0.,  1., 0.], 0.01)
fb9             = m.outer_faces_with_direction([0., 0., 0., 0., -1., 0.], 0.01)
fb10            = m.outer_faces_with_direction([0., 0., 0., 0., 0.,  1.], 0.01)
fb11            = m.outer_faces_with_direction([0., 0., 0., 0., 0., -1.], 0.01)

# mark it as boundary #4X
m.set_region(40, fb0)# Face the plasma is entering on
m.set_region(41, fb1)# Face the plasma is exiting
m.set_region(42, fb2)# Side face
m.set_region(43, fb3)# Side face
m.set_region(44, fb4)# Side face
m.set_region(45, fb5)# Side face
m.set_region(46, fb6)# Velocity extreme
m.set_region(47, fb7)# Velocity extreme
m.set_region(48, fb8)# Velocity extreme
m.set_region(49, fb9)# Velocity extreme
m.set_region(50, fb10)# Velocity extreme
m.set_region(51, fb11)# Velocity extreme

Do you think this would be the problem? I saw something similar in some of the examples, where I got this.

Sincerely,
Eric Comstock

102

this is the C++ code for outer_faces_with_direction()

static void
outer_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out,
            const std::string &condition="")
{
  mesh_faces_by_pts_list lst;
  dal::bit_vector cvlst, checked_pids, rejected_pids;

  darray normal_vector = in.pop().to_darray();
  scalar_type angle = in.pop().to_scalar();
  scalar_type threshold = gmm::vect_norm2(normal_vector)*cos(angle);
  bgeot::base_node un;

  dim_type elm_dim = m.dim();
  if (in.remaining() && in.front().is_integer())
   elm_dim = dim_type(in.pop().to_integer());

  if (in.remaining()) cvlst = in.pop().to_bit_vector(&m.convex_index());
  else cvlst = m.convex_index();

  for (dal::bv_visitor ic(cvlst); !ic.finished(); ++ic) {
    if (m.structure_of_convex(ic)->dim() == elm_dim) {
      for (short_type f = 0; f < m.structure_of_convex(ic)->nb_faces(); f++) {
        bgeot::mesh_structure::ind_pt_face_ct pt
          = m.ind_points_of_face_of_convex(ic, f);
        std::vector<size_type> p(pt.begin(), pt.end());
        size_type idx = lst.add_norepeat(mesh_faces_by_pts_list_elt(ic,f,p));
        lst[idx].cnt++;
      }
    }
  }
  size_type fcnt = 0;
  for (size_type i = 0; i < lst.size(); i++) {
    if (lst[i].cnt == 1) {
      un = m.mean_normal_of_face_of_convex(lst[i].cv, lst[i].f);
      if (gmm::vect_sp(normal_vector, un) < threshold - 1E-8)
        lst[i].cnt = -1;
      if (lst[i].cnt == 1) fcnt++;
    }
  }

  iarray w = out.pop().create_iarray(2,unsigned(fcnt));
  fcnt = 0;
  for (size_type j=0; j < lst.size(); j++) {
    if (lst[j].cnt == 1) {
      w(0,fcnt) = int(lst[j].cv+config::base_index());
      w(1,fcnt) = int(short_type(lst[j].f+config::base_index()));
      fcnt++;
    }
  }

}

check every single of fb0, fb1, etc for correctness.

Here you can find all mesh object methods that you can use for introspection

1 Like
103

Hello Kostas,

Thank you! I was not entirely sure what was going on in the code you sent me, so I commented it and put my best guesses as to what I thought was happening.

Do you know if any sort of list of getFEM-defined classes exists anywhere? The classes and their dependencies were very hard for me to parse.


static void                                                                    
outer_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out,           \\ Function that adds a region to the getfem::mesh object
            const std::string &condition="")                                   
{                                                                              
  mesh_faces_by_pts_list lst;                                                  \\ Define lst as a dal::dynamic_tree_sorted, which has the alias mesh_faces_by_pts_list. I could not find dynamic_tree_sorted in the documentation
  dal::bit_vector cvlst, checked_pids, rejected_pids;                          \\ Each of these is a bit_vector, which represents all the indexes of valid points of a mesh.
                                                                               
  darray normal_vector = in.pop().to_darray();                                 \\ class darray from getfemint.h; contains a single garray, which i am not sure what that is - could not be found in documentation. Is this an actual normal vector, or something else? Is this what is input in Python?
  scalar_type angle = in.pop().to_scalar();                                    \\ scalar_type is definitely a scalar. Is it what you set using --with-float and similar commands on the ./configure? I am not using that now, but I explored it last month as a potential speedup avenue. angle seems to be removing the top of the mexargs_in stack of arguments, which is presumably an angle.
  scalar_type threshold = gmm::vect_norm2(normal_vector)*cos(angle);           \\ Why is the norm of a normal vector being taken? Is that not by definition 1?
  bgeot::base_node un;                                                         \\ This seems to be a vertex of one of the mesh elements. It is presumably uninitialized.
                                                                               
  dim_type elm_dim = m.dim();                                                  \\ elm_dim is number of dimensions. dim_type is presumably a class or struct with a single public integer in it
  if (in.remaining() && in.front().is_integer())                               \\ Runs if there are remaining arguments in the stack, and the next one is an integer.
   elm_dim = dim_type(in.pop().to_integer());                                  \\ Why does this change the dimension?
                                                                               
  if (in.remaining()) cvlst = in.pop().to_bit_vector(&m.convex_index());       \\ If any arguments are left, the cvlst uninitialized list of indexes of valid points of a mesh is initialized with the next argument.
  else cvlst = m.convex_index();                                               \\ If not, all convexes in the mesh are initialized to cvlst
                                                                               
  for (dal::bv_visitor ic(cvlst); !ic.finished(); ++ic) {                      \\ loops through a "dal::bv_visitor" (which I could not find documentation for). Presumably it starts fully initialized to something related to cvlst, and then iterates to its completion - presumably whever ic.finished() is true.
    if (m.structure_of_convex(ic)->dim() == elm_dim) {                         \\ Checks to see if the current dimension (as measured by running the dim() function on whatever is at the memory location of m.structure_of_convex(ic)) is equal to the elm_dim variable, which I think is the total dimensions. I would geuss that m.structure_of_convex(ic) measures how many dimensions cvlst has.
      for (short_type f = 0; f < m.structure_of_convex(ic)->nb_faces(); f++) { \\ If dimensions match, start another for loop, which ends when the iterated variable f becomes greater than or equal to the number of faces of the convex (would that be 12 for my 6D case, or the number of 5-simplices on my boundary (probably around 100,000)?). f is a face id.
        bgeot::mesh_structure::ind_pt_face_ct pt                               \\ pt is really a gmm::tab_ref_index_ref, which is an indexed array (I think)
          = m.ind_points_of_face_of_convex(ic, f);                             \\ and it is initialized to the global point numbers for the convex ic.
        std::vector<size_type> p(pt.begin(), pt.end());                        \\ This returns a pointer of some kind, but I am not sure where it comes from or is going.
        size_type idx = lst.add_norepeat(mesh_faces_by_pts_list_elt(ic,f,p));  \\ Add a point to lst
        lst[idx].cnt++;                                                        \\ Increment lst's id list?
      }                                                                        
    }                                                                          
  }                                                                            \\ I think this for loop is supposed to populate lst with points
  size_type fcnt = 0;                                                          \\
  for (size_type i = 0; i < lst.size(); i++) {                                 \\ for loop iterated from 0 to lst.size(), which is presumably the total number of points within lst
    if (lst[i].cnt == 1) {                                                     \\ check if the point is in lst?
      un = m.mean_normal_of_face_of_convex(lst[i].cv, lst[i].f);               \\ Is this normal vector a normal vector?
      if (gmm::vect_sp(normal_vector, un) < threshold - 1E-8)                  \\ Check to see if it is closer to aligned than not
        lst[i].cnt = -1;                                                       \\ Remove from lst???
      if (lst[i].cnt == 1) fcnt++;                                             \\ increment fcnt - presumablt the number of points on the face
    }                                                                          
  }                                                                            
                                                                               
  iarray w = out.pop().create_iarray(2,unsigned(fcnt));                        \\ Is this also just a class with a single garray?
  fcnt = 0;                                                                    \\
  for (size_type j=0; j < lst.size(); j++) {                                   \\ for loop iterated from 0 to lst.size(), which is presumably the total number of points within lst
    if (lst[j].cnt == 1) {                                                     \\ check if the point is in lst?
      w(0,fcnt) = int(lst[j].cv+config::base_index());                         \\ If so, then the w array 
      w(1,fcnt) = int(short_type(lst[j].f+config::base_index()));              \\
      fcnt++;                                                                  \\ Increment fnct by 1
    }                                                                          
  }                                                                            
                                                                               
}

Where is m modified?

Sincerely,
Eric Comstock

104

it is not the function adding a region to a mesh, it is the function which is called when you do this:

it returns the array w which is what you see in fb0.

Not entirely well maintained but we have doxygen documentation
https://fossies.org/dox/getfem-5.4.4/namespacepython_1_1getfem.html

105

Hello Kostas,

Wait - if it generates the fb0 array, then why is the function a void? What is the output? Where does it come from?

Thank you for the class documentation! It is very helpful!

I did some testing on my region array sizes, and I got the following:

2025-06-17 04:33:43,416 [DEBUG] fb0 size: 61479
2025-06-17 04:33:43,417 [DEBUG] fb1 size: 48357
2025-06-17 04:33:43,417 [DEBUG] fb2 size: 60507
2025-06-17 04:33:43,418 [DEBUG] fb3 size: 47385
2025-06-17 04:33:43,418 [DEBUG] fb4 size: 61479
2025-06-17 04:33:43,419 [DEBUG] fb5 size: 48357
2025-06-17 04:33:43,420 [DEBUG] fb6 size: 61722
2025-06-17 04:33:43,420 [DEBUG] fb7 size: 47871
2025-06-17 04:33:43,421 [DEBUG] fb8 size: 62451
2025-06-17 04:33:43,421 [DEBUG] fb9 size: 47871
2025-06-17 04:33:43,422 [DEBUG] fb10 size: 18477
2025-06-17 04:33:43,422 [DEBUG] fb11 size: 18459

Does this count nodes multiple times for each time they appear in an element, or does it count elements, or what does it measure? I would expect these to all be the same.

Sincerely,
Eric Comstock

106

it is complicated :slightly_smiling_face:

It is element faces in your case. In your 6D mesh, fb0 is an array of 61479 5D faces. If you believe this is wrong, then we should debug the outer_faces_with_direction function.

The first number in each fb0 entry is the convex (i.e. element) number, and the second number is the local face number in the element.

107

Hello Kostas,

Okay - unfortunately, there is no good way for me to check that on a simplex mesh.

I will see about a temporary switch to a hexarractal mesh again (6D hyperrectangular), since for that it should be 1024 per face. Will this function work the same way in this case?

Sincerely,
Eric Comstock

108

Hello Kostas,

Okay - it seems that the regions and multipliers are the correct size when a hexarractal mesh is used. I assume they were way larger with a simplex mesh because many more simplices fit inside a hexarract.

The problem is that the exact same delay and error (INCTL(14)) are still appearing.

Sincerely,
Eric Comstock

109

ICNTL error is the symptom not the root cause.

The root cause is the dense tangent matrix.

The null-space elimination that the add_multiplier method adds to the model, solves the root cause but in a very costly way.

After having checked the region definition, you can check if the number of dofs in the region is also correct.

You can do this with something like print(mf.basic_dof_on_region(40).shape) for the entering plasma side.

you can try with different definitions of mf

mf = gf.MeshFem(m,1) #scalar field
mf.set_classical_fem(1) # 6-linear

or

mf.set_classical_fem(0) # discontinuous element-wise constant
110

Hello Kostas,

Okay - thank you! I tried what you suggested, and got a really confusing result.

2025-06-17 09:24:20,392 [DEBUG] fb0 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb1 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb2 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb3 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb4 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb5 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb6 size: 243
2025-06-17 09:24:20,393 [DEBUG] fb7 size: 243
2025-06-17 09:24:20,394 [DEBUG] fb8 size: 243
2025-06-17 09:24:20,394 [DEBUG] fb9 size: 243
2025-06-17 09:24:20,394 [DEBUG] fb10 size: 243
2025-06-17 09:24:20,394 [DEBUG] fb11 size: 243
2025-06-17 09:24:20,599 [DEBUG] fb0 shape: (0,)
2025-06-17 09:24:20,600 [DEBUG] fb1 shape: (0,)
2025-06-17 09:24:20,600 [DEBUG] fb2 shape: (0,)
2025-06-17 09:24:20,601 [DEBUG] fb3 shape: (0,)
2025-06-17 09:24:20,602 [DEBUG] fb4 shape: (0,)
2025-06-17 09:24:20,602 [DEBUG] fb5 shape: (0,)
2025-06-17 09:24:20,602 [DEBUG] fb6 shape: (0,)
2025-06-17 09:24:20,603 [DEBUG] fb7 shape: (0,)
2025-06-17 09:24:20,603 [DEBUG] fb8 shape: (0,)
2025-06-17 09:24:20,603 [DEBUG] fb9 shape: (0,)
2025-06-17 09:24:20,603 [DEBUG] fb10 shape: (0,)
2025-06-17 09:24:20,603 [DEBUG] fb11 shape: (0,)

For this case, the 243 makes sense (3^5, 4096 point test mesh has 4 points and 3 elements per side). I tried mf.basic_dof_on_region(40).shape() as well, and got the same result. Is mf.basic_dof_on_region(40) supposed to be a numpy array? I am not sure how something can have a size of 243 but a a shape of (0,).

I am working on trying different mf definitions currently.

Sincerely,
Eric Comstock

111

mf.basic_dof_on_region(40) and fb0 are two totally different things

112

Hello Kostas,

Right - apologies. I think I had my code in the wrong order.

I am now getting

2025-06-17 09:36:33,453 [DEBUG] fb0 size: 243
2025-06-17 09:36:33,454 [DEBUG] fb1 size: 243
2025-06-17 09:36:33,455 [DEBUG] fb2 size: 243
2025-06-17 09:36:33,456 [DEBUG] fb3 size: 243
2025-06-17 09:36:33,456 [DEBUG] fb4 size: 243
2025-06-17 09:36:33,457 [DEBUG] fb5 size: 243
2025-06-17 09:36:33,457 [DEBUG] fb6 size: 243
2025-06-17 09:36:33,457 [DEBUG] fb7 size: 243
2025-06-17 09:36:33,458 [DEBUG] fb8 size: 243
2025-06-17 09:36:33,458 [DEBUG] fb9 size: 243
2025-06-17 09:36:33,458 [DEBUG] fb10 size: 243
2025-06-17 09:36:33,459 [DEBUG] fb11 size: 243
2025-06-17 09:56:17,770 [DEBUG] region 0 shape: (1024,)
2025-06-17 09:56:17,771 [DEBUG] region 1 shape: (1024,)
2025-06-17 09:56:17,772 [DEBUG] region 2 shape: (1024,)
2025-06-17 09:56:17,773 [DEBUG] region 3 shape: (1024,)
2025-06-17 09:56:17,774 [DEBUG] region 4 shape: (1024,)
2025-06-17 09:56:17,775 [DEBUG] region 5 shape: (1024,)
2025-06-17 09:56:17,775 [DEBUG] region 6 shape: (1024,)
2025-06-17 09:56:17,776 [DEBUG] region 7 shape: (1024,)
2025-06-17 09:56:17,777 [DEBUG] region 8 shape: (1024,)
2025-06-17 09:56:17,778 [DEBUG] region 9 shape: (1024,)
2025-06-17 09:56:17,779 [DEBUG] region 10 shape: (1024,)
2025-06-17 09:56:17,780 [DEBUG] region 11 shape: (1024,)

1024 should be the number of vertices on my faces, and 243 should be my number of element faces on my face, so I think this is correct.

I will go back over the code for outer_faces, and see if I can glean anything.

The list of model variables and data is also different now that I am back to hexarracts:

message from gf_model_get follow:
List of model variables and data:
Variable       f                              1 copy   fem dependant     4096 doubles.
Variable       mult40                         1 copy   fem dependant     1024 doubles.
Variable       mult46                         1 copy   fem dependant     1024 doubles.
Variable       mult47                         1 copy   fem dependant     1024 doubles.
Variable       mult48                         1 copy   fem dependant     1024 doubles.
Variable       mult49                         1 copy   fem dependant     1024 doubles.
Variable       mult50                         1 copy   fem dependant     1024 doubles.
Variable       mult51                         1 copy   fem dependant     1024 doubles.
Data           B1                             1 copy   fem dependant     4096 doubles.
Data           B2                             1 copy   fem dependant     4096 doubles.
Data           B3                             1 copy   fem dependant     4096 doubles.
Data           DirichletData                  1 copy   fem dependant     4096 doubles.
Data           E1                             1 copy   fem dependant     4096 doubles.
Data           E2                             1 copy   fem dependant     4096 doubles.
Data           E3                             1 copy   fem dependant     4096 doubles.
Data           zeros                          1 copy   fem dependant     4096 doubles.

Sincerely,
Eric Comstock

113

when you work with hexarracts, how does the add_multiplier version compare to the add_filtered_fem_variable version of your program? do they produce the same numbers of dofs for all variables involved?

114

Hello Kostas,

Yes, I get

message from gf_model_get follow:
List of model variables and data:
Variable       f                              1 copy   fem dependant     4096 doubles.
Variable       mult40                         1 copy   fem dependant     1024 doubles.
Variable       mult46                         1 copy   fem dependant     1024 doubles.
Variable       mult47                         1 copy   fem dependant     1024 doubles.
Variable       mult48                         1 copy   fem dependant     1024 doubles.
Variable       mult49                         1 copy   fem dependant     1024 doubles.
Variable       mult50                         1 copy   fem dependant     1024 doubles.
Variable       mult51                         1 copy   fem dependant     1024 doubles.
Data           B1                             1 copy   fem dependant     4096 doubles.
Data           B2                             1 copy   fem dependant     4096 doubles.
Data           B3                             1 copy   fem dependant     4096 doubles.
Data           DirichletData                  1 copy   fem dependant     4096 doubles.
Data           E1                             1 copy   fem dependant     4096 doubles.
Data           E2                             1 copy   fem dependant     4096 doubles.
Data           E3                             1 copy   fem dependant     4096 doubles.
Data           zeros                          1 copy   fem dependant     4096 doubles.

and I still get the same ICNTL(14) error. That is why I originally switched from hexarracts to simplices.

I am doing some digging around in the reduction and extension matrices to try and see what the problem is.

Sincerely,
Eric Comstock

115

Hello Kostas,

Okay - so for both the hexarract and simplex cases, the total dofs (named V in Build a finite element method on a mesh — GetFEM) is a 4096 column vector.

In the hexarract case, the reduction and extension matrices are both binary matrices of rank 1024 (and can be rearranged into a matrix of

 I_1024  O
[         ]
 O       O

), and for the simplex cases, they are both of rank 3056 in mult46, which I tested, and can be rearranged into a matrix of

 I_3056  O
[         ],
 O       O

where I_N is the identity matrix of size N.

Both matrices seem to have a large null space, which is what selects the region in my cases. Is this supposed to happen beforehand (in V) or is it intended?

Sincerely,
Eric Comstock

116

let’s stick to hexarracts then.

for add_filtered_fem_variable, your model has 4096 + 7*1024 dofs

How many dofs does the model with add_multiplier have?

117

Hello Kostas,

I got exactly the same dofs, and the same List of model variables and data: results for both add_filtered_fem_variable and add_multiplier for the hexarract case.

Sincerely,
Eric Comstock

118

then both cases must be giving the ICNTL(14) error, right?

119

Hello Kostas,

Right - I get the ICNTL(14) when I use

  • Hexarracts with add_multiplier
  • Hexarracts with add_filtered_fem_variable
  • Simplices with add_filtered_fem_variable

I do not get the error when I use simplices with add_multiplier.

Sincerely,
Eric Comstock

120

let’s forget about the simplices, as if it does not exist

we only talk Hexarracts

run

import scipy.sparse.linalg
#import scipy.io

...

md.assembly("build_all")
K = md.tangent_matrix()
K = scipy.sparse.csc_matrix((K.csc_val(),*(K.csc_ind()[::-1])))
#scipy.io.savemat(f"K_{i}.mat", {'K':K}) # for exporting the matrix in matlab format
# K = scipy.io.loadmat(f"K_{i}.mat")["K"] # for importing the matrix from matlab format
print(K.getnnz())

...

and use matplotlib.pyplot.spy — Matplotlib 3.10.3 documentation in order to view the sparsity pattern of the matrix

1 Like