New Hex element

[daniele.barbera86]

Dear all,

recently we introduced the hex element, and I would like to know which modules currently support it, since I was planning to do my simulation of the flange using then. It looks like mesh is strongly affecting the results and hex could be the solution.

In addition there is any major change in the way we create our models in CUBIT? It could a nice way to start a small documentation on how to use them.

Thank you very much!

Daniele

[Karol Lewandowski]

There is no major change in Cubit, just use hex elements.

Currently contact tutorial supports them. They will be added into multifield in a couple of days as well. 

[Lukasz Kaczmraczyk]

Daniele,

Most of the codes should support Hex. You have to be careful only what approximation base you using, i.e. hex has all energetic spaces only for DEMKOWICZ_JACOBI_BASE. For example look at how is modified contact tutorial is, i.e. adv-0.

Depending on how code is implemented, sometimes is necessary to push operators for integration weights, when faces of hex are not parallel. That is also done in adv-0 example, look at functions, for volume

      auto jac_ptr = boost::make_shared<MatrixDouble>();
      auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
      auto det_ptr = boost::make_shared<VectorDouble>();
      pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
      pipeline.push_back(
          new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
      pipeline.push_back(
          new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
      pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
      pipeline.push_back(new OpSetHOInvJacToScalarBases(H1, inv_jac_ptr));
      pipeline.push_back(new OpSetHOWeights(det_ptr))

and for boundary,

      pipeline.push_back(new OpSetHOWeigthsOnFace())

​

You work in a cubit, as usual, only make hex mesh. It is a bit tricky, but simple geometry should work map. Ask questions if you run into problems. The best way is to try first.

Regards,

L

[Karol Lewandowski]

Lukasz, what should put I for example in cases when I search in multi index containers for MBTET? Should I search twice, for MBHEX as well? 

[Lukasz Kaczmraczyk]

Hi,

Depending on the container, but you can search for a range,

auto lo_it = some_multi_index_container.get<Ent_mi_tag>().lower_bound(get_id_for_min_type(MBTET ));

auto hi_it = some_multi_index_container.get<Ent_mi_tag>().lower_bound(get_id_for_min_type(MBHEX ))

for(;lo_it!=hi_it;lo_it) {

};

If you like to do it for arbitrary dimension,

const int dim = 3;

auto lo_it = some_multi_index_container.get<Ent_mi_tag>().lower_bound(get_id_for_min_type(moab::CN::TypeDimensionMap[dim].first));

auto hi_it = some_multi_index_container.get<Ent_mi_tag>().lower_bound(get_id_for_min_type(moab::CN::TypeDimensionMap[dim].second));

for(;lo_it!=hi_it;lo_it) {

};

L.

[Lukasz Kaczmraczyk]

small fix, should be,

auto hi_it = some_multi_index_container.get<Ent_mi_tag>().upper_bound(get_id_for_max_type(MBHEX ));

or 

auto hi_it = some_multi_index_container.get<Ent_mi_tag>().upper_bound(get_id_for_max_type(moab::CN::TypeDimensionMap[dim].second));

L.

[Karol Lewandowski]

Where is MESH_NODE_POSITIONS field? How do you post-process this? 

[Lukasz Kaczmraczyk]

Karol,

Where? In what context, lukasz/rotation_bc? 

If you asking bout that, I used the name "HO_POSITIONS" instead "MESH_NODE_POSITIONS" which is in fact not very precise. If you ask about contact, there is no HO geometry there, but you have to anyway use it for hex if the hex is not regular. In that case, jacobian is not constant, so integration weights are updated. 

L.

[Karol Lewandowski]

Ok, so you deleted HO from contact. Which updated module/tutorial is using it now?

[Lukasz Kaczmraczyk]

Then look to branch lukasz/rotation_bc and contact.

https://bitbucket.org/mofem/users-modules-cephas/src/c42085f52ca53991f58fee1afb1a355fde5bc49f/tutorials/adv-1/contact.cpp?at=lukasz%2Frotation_bc

[Karol Lewandowski]

Does it make sense to set 2nd order of "HO_POSITIONS" only on edges? 
How would you check if the geometry is 2nd order? 

CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);

CHKERR mField.get_moab().get_connectivity(&*tets.begin(), 1, nodes);

CHKERR mField.get_moab().get_number_entities_by_type(0, MBTET, nb_tets, true);

CHKERR mField.get_moab().get_number_entities_by_type(0, MBHEX, nb_hexes, true);

if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {

    data.is_ho_geometry = PETSC_TRUE;

}

[Lukasz Kaczmraczyk]

If you have TET, the second-order approximation is only on edges, that is why TET10, four nodes, and six edges. For HEX is a bit more complex, you can have elements that have nodes only on HEX edges, an element which has in addition noded in mid-quads faces of hax, and finally in the center of hex. Projection10NodeCoordsOnField hanldes now only edges. Need to be extended for other types. But it works OK with type from Cubit.

Projection10NodeCoordsOnField we have to as well to fix this for 2d.

[Karol Lewandowski]

So is better to just set 2nd order for all entities for now? For "HO_POSITIONS"? 

[Lukasz Kaczmraczyk]

In what context? Yes, we can transfer from Cubiut mesh up to 2d order, if this is what you mean. 

[Karol Lewandowski]

I mean this:

CHKERR simple->setFieldOrder("HO_POSITIONS", 2);

[Lukasz Kaczmraczyk]

Yes. Have more than 2nd makes sense if you calculate values for DOFs order higher than 2. Otherwise, those DOFs are only there, but the value is set to zero. Dead mass. Anyway, the projection will do only up to edges, but this is only information passed from Cubit. I do not if Cubit can handle orders bigger than 2. However, I think that we can extend that to centroids of Quads and Edges.

L.

[daniele.barbera86]

Thanks for all the replies. The mesh is ok since I worked always with hex and geometries are relatively simple. I am looking at adv-0 to run a simple case there.