How should we use bit level in MOFEM to manipulate prisms elements created during each load-step for contact mechanics module?

[zahur.ullah1]

How should we use bit level in MOFEM to manipulate prisms elements created during each load-step for contact mechanics module? 

I have already updated the mesh and have calculated new prisms elements, such as: 

// Set coordinate of nodes to current spatial procession, i.e. current configuration
ierr = contact_search_kd_tree.projectSptialPositionOnMesh("SPATIAL_POSITION"); CHKERRQ(ierr);
// Clear range of prisms elements
range_slave_master_prisms.clear();
// Create new prism for current configuration
ierr = contact_search_kd_tree.contactSearchAlgorithm(
 range_surf_master, range_surf_slave, contact_commondata_multi_index, range_slave_master_prisms
); CHKERRQ(ierr);

[Lukasz Kaczmraczyk]

Zahur,

Hi, Zahur,

First of all, you need to look how bit levels work, but you know that already, just for reference here is a link (Section Explaining BitRefLevel). 

So procure can work; like that:

1) you start with the case that all entities have bit-level set to 1000..0. That includes tetrahedra and prisms from the previous load step, let say n.

Tets:   00..01
Prisms: 00..01

2) next you set second bit to all tetrahedra:

Tets: 00..11

{
  Range tets_at_bit;
  BitRefLevel from_bit;
  from_bit.set(0); // 00..01
  ierr = m_mield.get_entities_by_type_and_ref_level(from_bit,BitRefLevel().set(),MBTET,tets_at_bit); CHKERRQ(ierr);
  from_bit.set(1); // 00..11
  ierr = m_field.seed_ref_level(tets_at_bit,from_bit,true);  CHKERRQ(ierr);
}

3) next you set last bit on all prisms, we will use last two bits as garbage,

Prisms: 00..01 set to 10..00

{
  Range prims_at_bit;
  BitRefLevel set_bit;

  set_bit.set(BITREFLEVEL_SIZE-1); // 10..00

  const RefEntity_multiIndex *ref_ent_ptr;

  ierr = m_field.get_ref_ents(&ref_ent_ptr);

  RefEntity_multiIndex::index<EntType_mi_tag>::type::iterator dit,hi_dit;

  dit = ref_ent_ptr->lower_bound(MBPRISM);

  hi_dit = ref_ent_ptr->upper_bound(MBPRISM);
  for(;dit!=hi_dit;dit++) {

    ref_ent_ptr->project<0>(dit),RefEntity_change_set_bit(set_bit)

  }

}

4) next you create new prisms, with your algorithm

// Set coordinate of nodes to current spatial procession, i.e. current configuration
ierr = contact_search_kd_tree.projectSptialPositionOnMesh("SPATIAL_POSITION"); CHKERRQ(ierr);
// Clear range of prisms elements
range_slave_master_prisms.clear();
// Create new prism for current configuration
ierr = contact_search_kd_tree.contactSearchAlgorithm(
 range_surf_master, range_surf_slave, contact_commondata_multi_index, range_slave_master_prisms
); CHKERRQ(ierr);

5)  Set second bit to new prisms

New prisms: 00..10

BitRefLevel set_bit;
set_bit.set(1); // 00..10
ierr = m_field. seed_ref_level(range_slave_master_prisms,set_bit,false); CHKERRQ(ierr);

5) Now you shift bit of all entities to the right as (this is to make sure we not run out of bit level as in MOFEM currently we have 32 bit level): 

Tets:       00..11 -> 00..01
New Prisms: 00..10 -> 00..01
Old Prims:  10..00 -> 01..00

ierr = mField.shift_right_bit_ref(1); CHKERRQ(ierr);

6) You can now delete all entities which have bit 01..00 (Old Prism). That should include fields & finite elements and entities in moab database. Do not do it now, first check if algorithm without deletion works. Note last two bits works as a garbage bin, you just accumulate on those bits everything with you do not need from previous load steps.

7) Now you have to create new DM and matrices with new adjacencies. Of course, if you build composite DM before, you only rebuild sub DM with contact part. Also, you can use block solver; field split solver. This is another step, you can do that later.

Now simply rebuild the whole problem and create vectors and matrices, 

ierr = DMDestroy(&dm);
ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
//set dm datastruture which created mofem datastructures
BitRefLevel problem_bit_level;
problem_bit_level.set(0); // 00..01
ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
//add elements to dm
ierr = DMMoFEMAddElement(dm,"ELASTIC"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
ierr = DMMoFEMAddElement(dm,"CONTACT_ELEM"); CHKERRQ(ierr);
ierr = DMSetUp(dm); CHKERRQ(ierr);

We have to clean interface with for setting bit levels, at the moment is a little bit messy and unintuitive. If I find time, I will do it.

Regards,

Lukasz

[zahur.ullah1]

How the bit-level enumerated in MOFEM, i.e. starting from left or right: So which one is correct: 

A:   10000...000

Or 

B:  00000...001

[Lukasz Kaczmraczyk]

B is correct. Can you fix my answer?

Should be like this,

1000..00 should be 00..0001
0100..00 should be 00..0010
1000..01 should be 00..0001

and the same for others.
 

[zahur.ullah1]

Hello Lukasz, 

After implementing your suggested procedure in my contact module. I found the following problem: 

New prisms (CONTACT_ELEM) are not included in the problem ELASTIC_MECHANICS. Currently, in the ELASTIC_MECHANICS we have only 513 elements which are the just summation of ELASTIC and FORCE_FE. 

name ELASTIC_MECHANICS Nb. elems 513

where 

FORCE_FE Nb. FEs 8

ELASTIC Nb. FEs 505

CONTACT_ELEM Nb. FEs 82

So what you suggest, do we need to remove old and add new prism to CONTACT_ELEM ?

And 

Do we need to update the fields, e.g. LAGMULT ?

Or 

There is something wrong with setting the bit-level for new prims. 

[Lukasz Kaczmraczyk]

Zahur, 

Look at step 5), is mistake, we should set second bit, instead third is on, should be

BitRefLevel set_bit;
set_bit.set(1); // 00..10 <-- was set_bit.set(2)

ierr = m_field. seed_ref_level(range_slave_master_prisms,set_bit,false); CHKERRQ(ierr);

See as well I changed step 3), simply sets bit to all existing prisms using low level functionality. This example shows that specialised interface with clear names is needed for solve problems like that.