255 {
256
259
260 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
261 if (fid < 0) {
263 "Unable to open file '%s'", file.c_str());
264 }
266 MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
267 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
268 if (verb > 1)
270 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
271 vf[1], vf[2],
v[0],
v[1],
v[2]);
272
273 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
275 "Cannot read MED file older than V2.2");
276 }
277
278 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
279 bzero(mesh_name, MED_NAME_SIZE + 1);
280 bzero(mesh_desc, MED_COMMENT_SIZE + 1);
281 med_int space_dim;
282 med_mesh_type mesh_type;
283 med_int mesh_dim, n_step;
284 char dt_unit[MED_SNAME_SIZE + 1];
285 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
286 med_sorting_type sorting_type;
287 med_axis_type axis_type;
288 if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
289 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
290 axis_name, axis_unit) < 0) {
292 "Unable to read mesh information");
293 }
294 if (verb > 0)
295 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Reading mesh %s nsteps %d", mesh_name,
296 n_step);
297
298 switch (axis_type) {
299 case MED_CARTESIAN:
300 break;
301 case MED_SPHERICAL:
302 case MED_CYLINDRICAL:
303 default:
305 "Curvilinear coordinate system implemented");
306 }
307 if (space_dim < 2) {
309 "Not implemented for space dim %d", space_dim);
310 }
311
313 {
314 MeshsetsManager *meshsets_manager_ptr;
315 CHKERR m_field.getInterface(meshsets_manager_ptr);
316 int max_id = 0;
318 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
319 }
320 max_id++;
322 std::string(mesh_name));
323 CubitMeshSet_multiIndex::index<
324 Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
325 cit =
326 meshsets_manager_ptr->getMeshsetsMultindex()
327 .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
329 max_id++;
330 mesh_meshset = cit->getMeshset();
332 }
333
334 med_bool change_of_coord, geo_transform;
335 med_int num_nodes = MEDmeshnEntity(
336 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
337 MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
338 if (num_nodes < 0) {
340 "Could not read number of MED nodes");
341 }
342 if (num_nodes == 0) {
344 "No nodes in MED mesh");
345 }
346 if (verb > 0)
347 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Read number of nodes %d", num_nodes);
348
349 std::vector<med_float> coord_med(space_dim * num_nodes);
350 if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
351 MED_NO_INTERLACE, &coord_med[0]) < 0) {
353 "Could not read MED node coordinates");
354 }
355
356
357 ReadUtilIface *iface;
358 vector<double *> arrays_coord;
360 CHKERR m_field.get_moab().query_interface(iface);
361 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
362 Range verts(startv, startv + num_nodes - 1);
363 std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
364 arrays_coord[0]);
365 std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
366 arrays_coord[1]);
367 if (space_dim == 2) {
368 std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
369 } else {
370 std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
371 arrays_coord[2]);
372 }
373 CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
374 family_elem_map.clear();
375
376
377 {
378 std::vector<med_int> nodes_tags(num_nodes, 0);
379 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
380 MED_NODE, MED_NO_GEOTYPE,
381 &nodes_tags[0]) < 0) {
382 nodes_tags.clear();
383 }
384 for (
int i = 0;
i < num_nodes;
i++) {
385 family_elem_map[nodes_tags.empty() ?
i : nodes_tags[
i]].insert(verts[
i]);
386 }
387 }
388
389
390 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
391
393
394 for (auto type : types) {
395
396
397 med_bool change_of_coord;
398 med_bool geo_transform;
399 med_int num_ele = MEDmeshnEntity(
400 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type,
401 MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
402
403
404 if (num_ele <= 0)
405 continue;
406
407 int num_nod_per_ele =
type % 100;
408 if (verb > 0)
410 << "Reading elements " << num_ele << " of type "
411 << moab::CN::EntityTypeName(ent_type) << " number of nodes "
412 << num_nod_per_ele;
413
414 std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
415 if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
416 MED_CELL, type, MED_NODAL,
417 MED_FULL_INTERLACE, &conn_med[0]) < 0) {
419 "Could not read MED elements");
420 }
421
423
424 if (ent_type != MBVERTEX) {
427 CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
428 starte, conn_moab);
429 switch (ent_type) {
430
431 case MBTET: {
432 int ii = 0;
433 for (int ee = 0; ee != num_ele; ee++) {
434 std::vector<EntityHandle>
n(num_nod_per_ele);
435 for (int nn = 0; nn != num_nod_per_ele; nn++) {
436 n[nn] = verts[conn_med[ii + nn] - 1];
437 }
438
439 std::array<int, 10> nodes_map{
440
441 1, 0, 2, 3,
442
443 4, 6, 5, 8,
444 7, 9
445
446 };
447
448 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
449 conn_moab[ii] =
n[nodes_map[nn]];
450 }
451 }
452 }
453 case MBHEX: {
454 int ii = 0;
455 for (int ee = 0; ee != num_ele; ee++) {
456 std::vector<EntityHandle>
n(num_nod_per_ele);
457 for (int nn = 0; nn != num_nod_per_ele; nn++) {
458 n[nn] = verts[conn_med[ii + nn] - 1];
459 }
460
461 std::array<int, 20> nodes_map{
462 0, 3, 2, 1, 4, 7, 6, 5,
463 11, 10, 9, 8,
464 16, 19, 18, 17,
465 15, 14, 13, 12
466
467 };
468
469 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
470 conn_moab[ii] =
n[nodes_map[nn]];
471 }
472 }
473 } break;
474 default: {
475 int ii = 0;
476 for (int ee = 0; ee != num_ele; ee++) {
477 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
478 conn_moab[ii] = verts[conn_med[ii] - 1];
479 }
480 }
481 }
482 }
483 CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
484 conn_moab);
485 ents =
Range(starte, starte + num_ele - 1);
486 } else {
487
488 int ii = 0;
489 std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
490 for (int ee = 0; ee != num_ele; ++ee)
491 for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
492 conn_moab[ii] = verts[conn_med[ii] - 1];
493 ents.insert_list(conn_moab.begin(), conn_moab.end());
494 }
495
496
497 CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
498
499
500 {
501 std::vector<med_int> fam(num_ele, 0);
502 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
503 MED_CELL, type, &fam[0]) < 0) {
505 "No family number for elements: using 0 as default family number";
506 }
507 for (
int j = 0;
j < num_ele;
j++) {
508 family_elem_map[fam[
j]].insert(ents[
j]);
509 }
510 }
511 }
512 }
513
514 if (MEDfileClose(fid) < 0) {
516 "Unable to close file '%s'", file.c_str());
517 }
518
520}
static std::vector< med_geometrie_element > moab2med_element_type(const EntityType type)
std::vector< EntityHandle > meshMeshsets
meshset for each mesh