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