MFEM  v4.6.0
Finite element discretization library
fmsconvert.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_FMS
15 
16 #include "fmsconvert.hpp"
17 #include <climits>
18 
19 using std::endl;
20 
21 // #define DEBUG_FMS_MFEM 1
22 // #define DEBUG_MFEM_FMS 1
23 
24 namespace mfem
25 {
26 
27 static inline int
28 FmsBasisTypeToMfemBasis(FmsBasisType b)
29 {
30  int retval = -1;
31  switch (b)
32  {
33  case FMS_NODAL_GAUSS_OPEN:
35  break;
36  case FMS_NODAL_GAUSS_CLOSED:
38  break;
39  case FMS_POSITIVE:
40  retval = mfem::BasisType::Positive;;
41  break;
42  case FMS_NODAL_UNIFORM_OPEN:
44  break;
45  case FMS_NODAL_UNIFORM_CLOSED:
47  break;
48  case FMS_NODAL_CHEBYSHEV_OPEN:
49  case FMS_NODAL_CHEBYSHEV_CLOSED:
50  mfem::out <<
51  "FMS_NODAL_CHEBYSHEV_OPEN, FMS_NODAL_CHEBYSHEV_CLOSED need conversion to MFEM types."
52  << endl;
53  break;
54  }
55  return retval;
56 }
57 
58 // The following function is unused (for now), so it is commented out to
59 // suppress compilation warning.
60 #if 0
61 /**
62 @brief Get the order and layout of the field.
63 */
64 static int
65 FmsFieldGetOrderAndLayout(FmsField f, FmsInt *f_order, FmsLayoutType *f_layout)
66 {
67  int err = 0;
68  FmsFieldDescriptor fd;
69  FmsLayoutType layout;
70  FmsScalarType data_type;
71  const void *data = nullptr;
72  FmsInt order = 0;
73 
74  FmsFieldGet(f, &fd, NULL, &layout, &data_type,
75  &data);
76 
77  FmsFieldDescriptorType f_fd_type;
78  FmsFieldDescriptorGetType(fd, &f_fd_type);
79  if (f_fd_type != FMS_FIXED_ORDER)
80  {
81  err = 1;
82  }
83  else
84  {
85  FmsFieldType field_type;
86  FmsBasisType basis_type;
87  FmsFieldDescriptorGetFixedOrder(fd, &field_type,
88  &basis_type, &order);
89  }
90 
91  *f_order = order;
92  *f_layout = layout;
93 
94  return err;
95 }
96 #endif
97 
98 /**
99 @brief This function converts an FmsField to an MFEM GridFunction.
100 @note I took some of the Pumi example code from the mesh conversion function
101  that converted coordinates and am trying to make it more general.
102  Coordinates are just another field so it seems like a good starting
103  point. We still have to support a bunch of other function types, etc.
104 */
105 template <typename DataType>
106 int
107 FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh,
108  GridFunction &func, bool setFE)
109 {
110  int err = 0;
111 
112  // NOTE: transplanted from the FmsMeshToMesh function
113  // We should do this work once and save it.
114  //--------------------------------------------------
115  FmsInt dim, n_elem, space_dim;
116 
117  // Find the first component that has coordinates - that will be the new mfem
118  // mesh.
119  FmsInt num_comp;
120  FmsMeshGetNumComponents(fms_mesh, &num_comp);
121  FmsComponent main_comp = NULL;
122  FmsField coords = NULL;
123  for (FmsInt comp_id = 0; comp_id < num_comp; comp_id++)
124  {
125  FmsComponent comp;
126  FmsMeshGetComponent(fms_mesh, comp_id, &comp);
127  FmsComponentGetCoordinates(comp, &coords);
128  if (coords) { main_comp = comp; break; }
129  }
130  if (!main_comp) { return 1; }
131  FmsComponentGetDimension(main_comp, &dim);
132  FmsComponentGetNumEntities(main_comp, &n_elem);
133  FmsInt n_ents[FMS_NUM_ENTITY_TYPES];
134  FmsInt n_main_parts;
135  FmsComponentGetNumParts(main_comp, &n_main_parts);
136  for (FmsInt et = FMS_VERTEX; et < FMS_NUM_ENTITY_TYPES; et++)
137  {
138  n_ents[et] = 0;
139  for (FmsInt part_id = 0; part_id < n_main_parts; part_id++)
140  {
141  FmsInt num_ents;
142  FmsComponentGetPart(main_comp, part_id, (FmsEntityType)et, NULL, NULL,
143  NULL, NULL, &num_ents);
144  n_ents[et] += num_ents;
145  }
146  }
147  //--------------------------------------------------
148 
149  // Interrogate the field.
150  FmsFieldDescriptor f_fd;
151  FmsLayoutType f_layout;
152  FmsScalarType f_data_type;
153  const void *f_data;
154  FmsFieldGet(f, &f_fd, &space_dim, &f_layout, &f_data_type,
155  &f_data);
156  // FmsFieldGet(coords, NULL, &space_dim, NULL, NULL, NULL);
157 
158  FmsInt f_num_dofs;
159  FmsFieldDescriptorGetNumDofs(f_fd, &f_num_dofs);
160 
161  // Access FMS data through this typed pointer.
162  auto src_data = reinterpret_cast<const DataType *>(f_data);
163 
164  FmsFieldDescriptorType f_fd_type;
165  FmsFieldDescriptorGetType(f_fd, &f_fd_type);
166  if (f_fd_type != FMS_FIXED_ORDER)
167  {
168  return 9;
169  }
170  FmsFieldType f_field_type;
171  FmsBasisType f_basis_type;
172  FmsInt f_order;
173  FmsFieldDescriptorGetFixedOrder(f_fd, &f_field_type,
174  &f_basis_type, &f_order);
175 
176  if (f_field_type != FMS_CONTINUOUS && f_field_type != FMS_DISCONTINUOUS &&
177  f_field_type != FMS_HDIV)
178  {
179  return 10;
180  }
181 
182  int btype = FmsBasisTypeToMfemBasis(f_basis_type);
183  if (btype < 0)
184  {
185  mfem::out << "\tInvalid BasisType: " << BasisType::Name(btype) << std::endl;
186  return 11;
187  }
188 
189  //------------------------------------------------------------------
190  if (setFE)
191  {
192  // We could assemble a name based on fe_coll.hpp rules and pass to
193  // FiniteElementCollection::New()
194 
195  mfem::FiniteElementCollection *fec = nullptr;
196  switch (f_field_type)
197  {
198  case FMS_DISCONTINUOUS:
199  fec = new L2_FECollection(f_order, dim, btype);
200  break;
201  case FMS_CONTINUOUS:
202  fec = new H1_FECollection(f_order, dim, btype);
203  break;
204  case FMS_HDIV:
205  fec = new RT_FECollection(f_order, dim);
206  break;
207  case FMS_HCURL:
208  case FMS_DISCONTINUOUS_WEIGHTED:
209  MFEM_ABORT("Field types FMS_HCURL and FMS_DISCONTINUOUS_WEIGHTED"
210  " are not supported yet.");
211  break;
212  }
213 
214  int ordering = (f_layout == FMS_BY_VDIM) ? Ordering::byVDIM : Ordering::byNODES;
215  auto fes = new FiniteElementSpace(mesh, fec, space_dim, ordering);
216  func.SetSpace(fes);
217  }
218  //------------------------------------------------------------------
219  const FmsInt nstride = (f_layout == FMS_BY_VDIM) ? space_dim : 1;
220  const FmsInt vstride = (f_layout == FMS_BY_VDIM) ? 1 : f_num_dofs;
221 
222  // Data reordering to store the data into func.
223  if ((FmsInt)(func.Size()) != f_num_dofs*space_dim)
224  {
225  return 12;
226  }
227 
228  mfem::FiniteElementSpace *fes = func.FESpace();
229  const int vdim = fes->GetVDim();
230  const mfem::FiniteElementCollection *fec = fes->FEColl();
231  const int vert_dofs = fec->DofForGeometry(mfem::Geometry::POINT);
232  const int edge_dofs = fec->DofForGeometry(mfem::Geometry::SEGMENT);
233  const int tri_dofs = fec->DofForGeometry(mfem::Geometry::TRIANGLE);
234  const int quad_dofs = fec->DofForGeometry(mfem::Geometry::SQUARE);
235  const int tet_dofs = fec->DofForGeometry(mfem::Geometry::TETRAHEDRON);
236  const int hex_dofs = fec->DofForGeometry(mfem::Geometry::CUBE);
237  int ent_dofs[FMS_NUM_ENTITY_TYPES];
238  ent_dofs[FMS_VERTEX] = vert_dofs;
239  ent_dofs[FMS_EDGE] = edge_dofs;
240  ent_dofs[FMS_TRIANGLE] = tri_dofs;
241  ent_dofs[FMS_QUADRILATERAL] = quad_dofs;
242  ent_dofs[FMS_TETRAHEDRON] = tet_dofs;
243  ent_dofs[FMS_HEXAHEDRON] = hex_dofs;
244  FmsInt fms_dof_offset = 0;
245  int mfem_ent_cnt[4] = {0,0,0,0}; // mfem entity counters, by dimension
246  int mfem_last_vert_cnt = 0;
249  if (dim >= 2 && edge_dofs > 0)
250  {
251  mfem::Array<int> ev;
252  for (int i = 0; i < mesh->GetNEdges(); i++)
253  {
254  mesh->GetEdgeVertices(i, ev);
255  int id = mfem_edge.GetId(ev[0], ev[1]);
256  if (id != i) { return 13; }
257  }
258  }
259  if (dim >= 3 &&
260  ((n_ents[FMS_TRIANGLE] > 0 && tri_dofs > 0) ||
261  (n_ents[FMS_QUADRILATERAL] > 0 && quad_dofs > 0)))
262  {
263  mfem::Array<int> fv;
264  for (int i = 0; i < mesh->GetNFaces(); i++)
265  {
266  mesh->GetFaceVertices(i, fv);
267  if (fv.Size() == 3) { fv.Append(INT_MAX); }
268  // HashTable uses the smallest 3 of the 4 indices to hash Hashed4
269  int id = mfem_face.GetId(fv[0], fv[1], fv[2], fv[3]);
270  if (id != i) { return 14; }
271  }
272  }
273 
274  // Loop over all parts of the main component
275  for (FmsInt part_id = 0; part_id < n_main_parts; part_id++)
276  {
277  // Loop over all entity types in the part
278  for (FmsInt et = FMS_VERTEX; et < FMS_NUM_ENTITY_TYPES; et++)
279  {
280  FmsDomain domain;
281  FmsIntType ent_id_type;
282  const void *ents;
283  const FmsOrientation *ents_ori;
284  FmsInt num_ents;
285  FmsComponentGetPart(main_comp, part_id, (FmsEntityType)et, &domain,
286  &ent_id_type, &ents, &ents_ori, &num_ents);
287  if (num_ents == 0) { continue; }
288  if (ent_dofs[et] == 0)
289  {
290  if (et == FMS_VERTEX) { mfem_last_vert_cnt = mfem_ent_cnt[et]; }
291  mfem_ent_cnt[FmsEntityDim[et]] += num_ents;
292  continue;
293  }
294 
295  if (ents != NULL &&
296  (ent_id_type != FMS_INT32 && ent_id_type != FMS_UINT32))
297  {
298  return 15;
299  }
300  if (ents_ori != NULL)
301  {
302  return 16;
303  }
304 
305  if (et == FMS_VERTEX)
306  {
307  const int mfem_dof_offset = mfem_ent_cnt[0]*vert_dofs;
308  for (FmsInt i = 0; i < num_ents*vert_dofs; i++)
309  {
310  for (int j = 0; j < vdim; j++)
311  {
312  const int idx = i*nstride+j*vstride;
313  func(mfem_dof_offset*nstride+idx) =
314  static_cast<double>(src_data[fms_dof_offset*nstride+idx]);
315  }
316  }
317  fms_dof_offset += num_ents*vert_dofs;
318  mfem_last_vert_cnt = mfem_ent_cnt[et];
319  mfem_ent_cnt[0] += num_ents;
320  continue;
321  }
322  mfem::Array<int> dofs;
323  if (FmsEntityDim[et] == dim)
324  {
325  for (FmsInt e = 0; e < num_ents; e++)
326  {
327  fes->GetElementInteriorDofs(mfem_ent_cnt[dim]+e, dofs);
328  for (int i = 0; i < ent_dofs[et]; i++, fms_dof_offset++)
329  {
330  for (int j = 0; j < vdim; j++)
331  {
332  func(fes->DofToVDof(dofs[i],j)) =
333  static_cast<double>(src_data[fms_dof_offset*nstride+j*vstride]);
334  }
335  }
336  }
337  mfem_ent_cnt[dim] += num_ents;
338  continue;
339  }
340  const FmsInt nv = FmsEntityNumVerts[et];
341  mfem::Array<int> ents_verts(num_ents*nv), m_ev;
342  const int *ei = (const int *)ents;
343  if (ents == NULL)
344  {
345  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL, FMS_INT32,
346  0, ents_verts.GetData(), num_ents);
347  }
348  else
349  {
350  for (FmsInt i = 0; i < num_ents; i++)
351  {
352  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL,
353  FMS_INT32, ei[i], &ents_verts[i*nv], 1);
354  }
355  }
356  for (int i = 0; i < ents_verts.Size(); i++)
357  {
358  ents_verts[i] += mfem_last_vert_cnt;
359  }
360  const int *perm;
361  switch ((FmsEntityType)et)
362  {
363  case FMS_EDGE:
364  {
365  for (FmsInt part_ent_id = 0; part_ent_id < num_ents; part_ent_id++)
366  {
367  const int *ev = &ents_verts[2*part_ent_id];
368  int mfem_edge_id = mfem_edge.FindId(ev[0], ev[1]);
369  if (mfem_edge_id < 0)
370  {
371  return 17;
372  }
373  mesh->GetEdgeVertices(mfem_edge_id, m_ev);
374  int ori = (ev[0] == m_ev[0]) ? 0 : 1;
376  fes->GetEdgeInteriorDofs(mfem_edge_id, dofs);
377  for (int i = 0; i < edge_dofs; i++)
378  {
379  for (int j = 0; j < vdim; j++)
380  {
381  func(fes->DofToVDof(dofs[i],j)) =
382  static_cast<double>(src_data[(fms_dof_offset+perm[i])*nstride+j*vstride]);
383  }
384  }
385  fms_dof_offset += edge_dofs;
386  }
387  break;
388  }
389  case FMS_TRIANGLE:
390  {
391  for (FmsInt part_ent_id = 0; part_ent_id < num_ents; part_ent_id++)
392  {
393  const int *tv = &ents_verts[3*part_ent_id];
394  int mfem_face_id = mfem_face.FindId(tv[0], tv[1], tv[2], INT_MAX);
395  if (mfem_face_id < 0)
396  {
397  return 18;
398  }
399  mesh->GetFaceVertices(mfem_face_id, m_ev);
400  int ori = 0;
401  while (tv[ori] != m_ev[0]) { ori++; }
402  ori = (tv[(ori+1)%3] == m_ev[1]) ? 2*ori : 2*ori+1;
404  fes->GetFaceInteriorDofs(mfem_face_id, dofs);
405  for (int i = 0; i < tri_dofs; i++)
406  {
407  for (int j = 0; j < vdim; j++)
408  {
409  func(fes->DofToVDof(dofs[i],j)) =
410  static_cast<double>(src_data[(fms_dof_offset+perm[i])*nstride+j*vstride]);
411  }
412  }
413  fms_dof_offset += tri_dofs;
414  }
415  break;
416  }
417  case FMS_QUADRILATERAL:
418  {
419  for (FmsInt part_ent_id = 0; part_ent_id < num_ents; part_ent_id++)
420  {
421  const int *qv = &ents_verts[4*part_ent_id];
422  int mfem_face_id = mfem_face.FindId(qv[0], qv[1], qv[2], qv[3]);
423  if (mfem_face_id < 0) { return 19; }
424  mesh->GetFaceVertices(mfem_face_id, m_ev);
425  int ori = 0;
426  while (qv[ori] != m_ev[0]) { ori++; }
427  ori = (qv[(ori+1)%4] == m_ev[1]) ? 2*ori : 2*ori+1;
429  fes->GetFaceInteriorDofs(mfem_face_id, dofs);
430  for (int i = 0; i < quad_dofs; i++)
431  {
432  for (int j = 0; j < vdim; j++)
433  {
434  func(fes->DofToVDof(dofs[i],j)) =
435  static_cast<double>(src_data[(fms_dof_offset+perm[i])*nstride+j*vstride]);
436  }
437  }
438  fms_dof_offset += quad_dofs;
439  }
440  break;
441  }
442  default: break;
443  }
444  mfem_ent_cnt[FmsEntityDim[et]] += num_ents;
445  }
446  }
447 
448  return err;
449 }
450 
451 int
452 FmsMeshToMesh(FmsMesh fms_mesh, Mesh **mfem_mesh)
453 {
454  FmsInt dim, n_vert, n_elem, n_bdr_elem, space_dim;
455 
456  // Find the first component that has coordinates - that will be the new mfem
457  // mesh.
458  FmsInt num_comp;
459  FmsMeshGetNumComponents(fms_mesh, &num_comp);
460  FmsComponent main_comp = NULL;
461  FmsField coords = NULL;
462  for (FmsInt comp_id = 0; comp_id < num_comp; comp_id++)
463  {
464  FmsComponent comp;
465  FmsMeshGetComponent(fms_mesh, comp_id, &comp);
466  FmsComponentGetCoordinates(comp, &coords);
467  if (coords) { main_comp = comp; break; }
468  }
469  if (!main_comp) { return 1; }
470  FmsComponentGetDimension(main_comp, &dim);
471  FmsComponentGetNumEntities(main_comp, &n_elem);
472  FmsInt n_ents[FMS_NUM_ENTITY_TYPES];
473  FmsInt n_main_parts;
474  FmsComponentGetNumParts(main_comp, &n_main_parts);
475 
476 #define RENUMBER_ENTITIES
477 #ifdef RENUMBER_ENTITIES
478  // I noticed that to get domains working right, since they appear to be
479  // defined in a local vertex numbering scheme, we have to offset the vertex
480  // ids that MFEM makes for shapes to move them to the coordinates in the
481  // current domain.
482 
483  // However, parts would just be a set of element ids in the current domain
484  // and it does not seem appropriate to offset the points in that case.
485  // Should domains be treated specially?
486  int *verts_per_part = new int[n_main_parts];
487 #endif
488 
489  // Sum the counts for each entity type across parts.
490  for (FmsInt et = FMS_VERTEX; et < FMS_NUM_ENTITY_TYPES; et++)
491  {
492  n_ents[et] = 0;
493  for (FmsInt part_id = 0; part_id < n_main_parts; part_id++)
494  {
495  FmsInt num_ents;
496  FmsComponentGetPart(main_comp, part_id, (FmsEntityType)et, NULL, NULL,
497  NULL, NULL, &num_ents);
498  n_ents[et] += num_ents;
499 #ifdef RENUMBER_ENTITIES
500  if (et == FMS_VERTEX)
501  {
502  verts_per_part[part_id] = num_ents;
503  }
504 #endif
505  }
506  }
507  n_vert = n_ents[FMS_VERTEX];
508 
509 #ifdef RENUMBER_ENTITIES
510  int *verts_start = new int[n_main_parts];
511  verts_start[0] = 0;
512  for (FmsInt i = 1; i < n_main_parts; ++i)
513  {
514  verts_start[i] = verts_start[i-1] + verts_per_part[i-1];
515  }
516 #endif
517 
518  // The first related component of dimension dim-1 will be the boundary of the
519  // new mfem mesh.
520  FmsComponent bdr_comp = NULL;
521  FmsInt num_rel_comps;
522  const FmsInt *rel_comp_ids;
523  FmsComponentGetRelations(main_comp, &rel_comp_ids, &num_rel_comps);
524  for (FmsInt i = 0; i < num_rel_comps; i++)
525  {
526  FmsComponent comp;
527  FmsMeshGetComponent(fms_mesh, rel_comp_ids[i], &comp);
528  FmsInt comp_dim;
529  FmsComponentGetDimension(comp, &comp_dim);
530  if (comp_dim == dim-1) { bdr_comp = comp; break; }
531  }
532  if (bdr_comp)
533  {
534  FmsComponentGetNumEntities(bdr_comp, &n_bdr_elem);
535  }
536  else
537  {
538  n_bdr_elem = 0;
539  }
540 
541  FmsFieldGet(coords, NULL, &space_dim, NULL, NULL, NULL);
542  int err = 0;
543  Mesh *mesh = nullptr;
544  mesh = new Mesh(dim, n_vert, n_elem, n_bdr_elem, space_dim);
545 
546  // Element tags
547  FmsInt num_tags;
548  FmsMeshGetNumTags(fms_mesh, &num_tags);
549  FmsTag elem_tag = NULL, bdr_tag = NULL;
550  for (FmsInt tag_id = 0; tag_id < num_tags; tag_id++)
551  {
552  FmsTag tag;
553  FmsMeshGetTag(fms_mesh, tag_id, &tag);
554  FmsComponent comp;
555  FmsTagGetComponent(tag, &comp);
556  if (!elem_tag && comp == main_comp)
557  {
558 #if DEBUG_FMS_MFEM
559  const char *tn = NULL;
560  FmsTagGetName(tag, &tn);
561  mfem::out << "Found element tag " << tn << std::endl;
562 #endif
563  elem_tag = tag;
564  }
565  else if (!bdr_tag && comp == bdr_comp)
566  {
567 #if DEBUG_FMS_MFEM
568  const char *tn = NULL;
569  FmsTagGetName(tag, &tn);
570  mfem::out << "Found boundary tag " << tn << std::endl;
571 #endif
572  bdr_tag = tag;
573  }
574  }
575  FmsIntType attr_type;
576  const void *v_attr, *v_bdr_attr;
577  mfem::Array<int> attr, bdr_attr;
578  FmsInt num_attr;
579  // Element attributes
580  if (elem_tag)
581  {
582  FmsTagGet(elem_tag, &attr_type, &v_attr, &num_attr);
583  if (attr_type == FMS_UINT8)
584  {
585  mfem::Array<uint8_t> at((uint8_t*)v_attr, num_attr);
586  attr = at;
587  }
588  else if (attr_type == FMS_INT32 || attr_type == FMS_UINT32)
589  {
590  attr.MakeRef((int*)v_attr, num_attr);
591  }
592  else
593  {
594  err = 1; // "attribute type not supported!"
595  goto func_exit;
596  }
597  }
598  // Boundary attributes
599  if (bdr_tag)
600  {
601  FmsTagGet(bdr_tag, &attr_type, &v_bdr_attr, &num_attr);
602  if (attr_type == FMS_UINT8)
603  {
604  mfem::Array<uint8_t> at((uint8_t*)v_bdr_attr, num_attr);
605  bdr_attr = at;
606  }
607  else if (attr_type == FMS_INT32 || attr_type == FMS_UINT32)
608  {
609  bdr_attr.MakeRef((int*)v_bdr_attr, num_attr);
610  }
611  else
612  {
613  err = 2; // "bdr attribute type not supported!"
614  goto func_exit;
615  }
616  }
617 
618  // Add elements
619  for (FmsInt part_id = 0; part_id < n_main_parts; part_id++)
620  {
621  for (int et = FMS_VERTEX; et < FMS_NUM_ENTITY_TYPES; et++)
622  {
623  if (FmsEntityDim[et] != dim) { continue; }
624 
625  FmsDomain domain;
626  FmsIntType elem_id_type;
627  const void *elem_ids;
628  const FmsOrientation *elem_ori;
629  FmsInt num_elems;
630  FmsComponentGetPart(main_comp, part_id, (FmsEntityType)et, &domain,
631  &elem_id_type, &elem_ids, &elem_ori, &num_elems);
632 
633  if (num_elems == 0) { continue; }
634 
635  if (elem_ids != NULL &&
636  (elem_id_type != FMS_INT32 && elem_id_type != FMS_UINT32))
637  {
638  err = 3; goto func_exit;
639  }
640  if (elem_ori != NULL)
641  {
642  err = 4; goto func_exit;
643  }
644 
645  const FmsInt nv = FmsEntityNumVerts[et];
646  mfem::Array<int> ents_verts(num_elems*nv);
647  if (elem_ids == NULL)
648  {
649  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL, FMS_INT32,
650  0, ents_verts.GetData(), num_elems);
651  }
652  else
653  {
654  const int *ei = (const int *)elem_ids;
655  for (FmsInt i = 0; i < num_elems; i++)
656  {
657  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL, FMS_INT32,
658  ei[i], &ents_verts[i*nv], 1);
659  }
660  }
661  const int elem_offset = mesh->GetNE();
662  switch ((FmsEntityType)et)
663  {
664  case FMS_EDGE:
665  err = 5;
666  goto func_exit;
667  break;
668  case FMS_TRIANGLE:
669 #ifdef RENUMBER_ENTITIES
670  // The domain vertices/edges were defined in local ordering. We
671  // now have a set of triangle vertices defined in terms of local
672  // vertex numbers. Renumber them to a global numbering.
673  for (FmsInt i = 0; i < num_elems*3; i++)
674  {
675  ents_verts[i] += verts_start[part_id];
676  }
677 #endif
678 
679  for (FmsInt i = 0; i < num_elems; i++)
680  {
681  mesh->AddTriangle(
682  &ents_verts[3*i], elem_tag ? attr[elem_offset+i] : 1);
683  }
684  break;
685  case FMS_QUADRILATERAL:
686 #ifdef RENUMBER_ENTITIES
687  for (FmsInt i = 0; i < num_elems*4; i++)
688  {
689  ents_verts[i] += verts_start[part_id];
690  }
691 #endif
692  for (FmsInt i = 0; i < num_elems; i++)
693  {
694  mesh->AddQuad(&ents_verts[4*i], elem_tag ? attr[elem_offset+i] : 1);
695  }
696  break;
697  case FMS_TETRAHEDRON:
698 #ifdef RENUMBER_ENTITIES
699  for (FmsInt i = 0; i < num_elems*4; i++)
700  {
701  ents_verts[i] += verts_start[part_id];
702  }
703 #endif
704  for (FmsInt i = 0; i < num_elems; i++)
705  {
706  mesh->AddTet(&ents_verts[4*i], elem_tag ? attr[elem_offset+i] : 1);
707  }
708  break;
709 
710  // TODO: What about wedges and pyramids?
711 
712 
713  case FMS_HEXAHEDRON:
714 #ifdef RENUMBER_ENTITIES
715  for (FmsInt i = 0; i < num_elems*8; i++)
716  {
717  ents_verts[i] += verts_start[part_id];
718  }
719 
720 #endif
721  for (FmsInt i = 0; i < num_elems; i++)
722  {
723  const int *hex_verts = &ents_verts[8*i];
724 #if 0
725  const int reorder[8] = {0, 1, 2, 3, 5, 4, 6, 7};
726  const int new_verts[8] = {hex_verts[reorder[0]],
727  hex_verts[reorder[1]],
728  hex_verts[reorder[2]],
729  hex_verts[reorder[3]],
730  hex_verts[reorder[4]],
731  hex_verts[reorder[5]],
732  hex_verts[reorder[6]],
733  hex_verts[reorder[7]]
734  };
735  hex_verts = new_verts;
736 #endif
737  mesh->AddHex(hex_verts, elem_tag ? attr[elem_offset+i] : 1);
738  }
739  break;
740  default:
741  break;
742  }
743  }
744  }
745 
746  // Add boundary elements
747  if (bdr_comp && n_bdr_elem > 0)
748  {
749  FmsInt n_bdr_parts;
750  FmsComponentGetNumParts(bdr_comp, &n_bdr_parts);
751 
752  for (FmsInt part_id = 0; part_id < n_bdr_parts; part_id++)
753  {
754  for (int et = FMS_VERTEX; et < FMS_NUM_ENTITY_TYPES; et++)
755  {
756  if (FmsEntityDim[et] != dim-1) { continue; }
757 
758  FmsDomain domain;
759  FmsIntType elem_id_type;
760  const void *elem_ids;
761  const FmsOrientation *elem_ori;
762  FmsInt num_elems;
763  FmsComponentGetPart(bdr_comp, part_id, (FmsEntityType)et, &domain,
764  &elem_id_type, &elem_ids, &elem_ori, &num_elems);
765  if (num_elems == 0) { continue; }
766 
767  if (elem_ids != NULL &&
768  (elem_id_type != FMS_INT32 && elem_id_type != FMS_UINT32))
769  {
770  err = 6; goto func_exit;
771  }
772  if (elem_ori != NULL)
773  {
774  err = 7; goto func_exit;
775  }
776 
777  const FmsInt nv = FmsEntityNumVerts[et];
778  mfem::Array<int> ents_verts(num_elems*nv);
779  if (elem_ids == NULL)
780  {
781  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL, FMS_INT32,
782  0, ents_verts.GetData(), num_elems);
783  }
784  else
785  {
786  const int *ei = (const int *)elem_ids;
787  for (FmsInt i = 0; i < num_elems; i++)
788  {
789  FmsDomainGetEntitiesVerts(domain, (FmsEntityType)et, NULL,
790  FMS_INT32, ei[i], &ents_verts[i*nv], 1);
791  }
792  }
793  const int elem_offset = mesh->GetNBE();
794  switch ((FmsEntityType)et)
795  {
796  case FMS_EDGE:
797  for (FmsInt i = 0; i < num_elems; i++)
798  {
799  mesh->AddBdrSegment(
800  &ents_verts[2*i], bdr_tag ? bdr_attr[elem_offset+i] : 1);
801  }
802  break;
803  case FMS_TRIANGLE:
804  for (FmsInt i = 0; i < num_elems; i++)
805  {
806  mesh->AddBdrTriangle(
807  &ents_verts[3*i], bdr_tag ? bdr_attr[elem_offset+i] : 1);
808  }
809  break;
810  case FMS_QUADRILATERAL:
811  for (FmsInt i = 0; i < num_elems; i++)
812  {
813  mesh->AddBdrQuad(
814  &ents_verts[4*i], bdr_tag ? bdr_attr[elem_offset+i] : 1);
815  }
816  break;
817  default:
818  break;
819  }
820  }
821  }
822  }
823 
824 #ifdef RENUMBER_ENTITIES
825  delete [] verts_per_part;
826  delete [] verts_start;
827 #endif
828 
829  // Transfer coordinates
830  {
831  // Set the vertex coordinates to zero
832  const double origin[3] = {0.,0.,0.};
833  for (FmsInt vi = 0; vi < n_vert; vi++)
834  {
835  mesh->AddVertex(origin);
836  }
837 
838  // Finalize the mesh topology
839  mesh->FinalizeTopology();
840 
841  FmsFieldDescriptor coords_fd = NULL;
842  FmsLayoutType coords_layout;
843  FmsFieldGet(coords, &coords_fd, NULL, &coords_layout, NULL, NULL);
844  if (!coords_fd)
845  {
846  mfem::err << "Error reading the FMS mesh coords' FieldDescriptor." << std::endl;
847  err = 8;
848  goto func_exit;
849  }
850  FmsInt coords_order = 0;
851  FmsBasisType coords_btype = FMS_NODAL_GAUSS_CLOSED;
852  FmsFieldType coords_ftype = FMS_CONTINUOUS;
853  FmsFieldDescriptorGetFixedOrder(coords_fd, &coords_ftype, &coords_btype,
854  &coords_order);
855  // Maybe this is extra but it seems mesh->SetCurvature assumes
856  // btype=1. Maybe protects us against corrupt data.
857  if (coords_btype != FMS_NODAL_GAUSS_CLOSED)
858  {
859  mfem::err << "Error reading FMS mesh coords." << std::endl;
860  err = 9;
861  goto func_exit;
862  }
863 
864  // Switch to mfem::Mesh with nodes (interpolates the linear coordinates)
865  const bool discont = (coords_ftype == FMS_DISCONTINUOUS);
866  mesh->SetCurvature(coords_order, discont, space_dim,
867  (coords_layout == FMS_BY_VDIM) ?
869 
870  // Finalize mesh construction
871  mesh->Finalize();
872 
873  // Set the high-order mesh nodes
874  mfem::GridFunction &nodes = *mesh->GetNodes();
875  FmsFieldToGridFunction<double>(fms_mesh, coords, mesh, nodes, false);
876  }
877 
878 func_exit:
879 
880  if (err)
881  {
882  delete mesh;
883  }
884  else
885  {
886  *mfem_mesh = mesh;
887  }
888  return err;
889 }
890 
891 bool
892 BasisTypeToFmsBasisType(int bt, FmsBasisType &btype)
893 {
894  bool retval = false;
895  switch (bt)
896  {
898  // mfem::out << "mfem::BasisType::GaussLegendre -> FMS_NODAL_GAUSS_OPEN" << endl;
899  btype = FMS_NODAL_GAUSS_OPEN;
900  retval = true;
901  break;
903  // mfem::out << "mfem::BasisType::GaussLobato -> FMS_NODAL_GAUSS_CLOSED" << endl;
904  btype = FMS_NODAL_GAUSS_CLOSED;
905  retval = true;
906  break;
908  // mfem::out << "mfem::BasisType::Positive -> FMS_POSITIVE" << endl;
909  btype = FMS_POSITIVE;
910  retval = true;
911  break;
913  // mfem::out << "mfem::BasisType::OpenUniform -> ?" << endl;
914  btype = FMS_NODAL_UNIFORM_OPEN;
915  retval = true;
916  break;
918  // mfem::out << "mfem::BasisType::ClosedUniform -> ?" << endl;
919  btype = FMS_NODAL_UNIFORM_CLOSED;
920  retval = true;
921  break;
923  // mfem::out << "mfem::BasisType::OpenHalfUniform -> ?" << endl;
924  break;
926  // mfem::out << "mfem::BasisType::Serendipity -> ?" << endl;
927  break;
929  // mfem::out << "mfem::BasisType::ClosedGL -> ?" << endl;
930  break;
931 
932  }
933  /*
934  Which MFEM types map to:?
935  FMS_NODAL_CHEBYSHEV_OPEN,
936  FMS_NODAL_CHEBYSHEV_CLOSED,
937  */
938 
939  return retval;
940 }
941 
942 /**
943 @note We add the FMS field descriptor and field in here so we can only do it
944  after successfully validating the inputs (handling certain grid function
945  types, etc.)
946 */
947 int
948 GridFunctionToFmsField(FmsDataCollection dc,
949  FmsComponent comp,
950  const std::string &fd_name,
951  const std::string &field_name,
952  const Mesh *mesh,
953  const GridFunction *gf,
954  FmsField *outfield)
955 {
956  if (!dc) { return 1; }
957  if (!comp) { return 2; }
958  if (!mesh) { return 3; }
959  if (!gf) { return 4; }
960  if (!outfield) { return 5; }
961 
962  double *c = gf->GetData();
963 
964  const mfem::FiniteElementSpace *fespace = gf->FESpace();
965  const mfem::FiniteElementCollection *fecoll = fespace->FEColl();
966 
967 #ifdef DEBUG_MFEM_FMS
968  mfem::out << "Adding FMS field for " << field_name << "..." << endl;
969 #endif
970 
971  /* Q: No getter for the basis, do different kinds of FECollection have
972  implied basis? There are two subclasses that actually have the getter,
973  maybe those aren't implied? */
974  FmsInt order = 1;
975  int vdim = 1;
976  FmsFieldType ftype = FMS_CONTINUOUS;
977  FmsBasisType btype = FMS_NODAL_GAUSS_CLOSED;
978  switch (fecoll->GetContType())
979  {
981  {
982  ftype = FMS_CONTINUOUS;
983  order = static_cast<FmsInt>(fespace->GetOrder(0));
984  vdim = gf->VectorDim();
985  auto fec = dynamic_cast<const mfem::H1_FECollection *>(fecoll);
986  if (fec != nullptr)
987  {
988  if (!BasisTypeToFmsBasisType(fec->GetBasisType(), btype))
989  {
990  mfem::err << "Error converting MFEM basis type to FMS for"
991  " FMS_CONTINUOUS." << std::endl;
992  return 6;
993  }
994  }
995  break;
996  }
998  {
999  ftype = FMS_DISCONTINUOUS;
1000  order = static_cast<FmsInt>(fespace->GetOrder(0));
1001  vdim = gf->VectorDim();
1002  auto fec = dynamic_cast<const mfem::L2_FECollection *>(fecoll);
1003  if (fec != nullptr)
1004  {
1005  if (!BasisTypeToFmsBasisType(fec->GetBasisType(), btype))
1006  {
1007  mfem::err << "Error converting MFEM basis type to FMS for"
1008  " FMS_DISCONTINUOUS." << std::endl;
1009  return 7;
1010  }
1011  }
1012  break;
1013  }
1015  {
1016  mfem::out << "Warning, unsupported ContType (TANGENTIAL) for "
1017  << field_name << ". Using FMS_CONTINUOUS." << std::endl;
1018  break;
1019  }
1021  {
1022  ftype = FMS_HDIV;
1023  // This RT_FECollection type seems to arise from "RT" fields such as "RT_3D_P1".
1024  // Checking fe_coll.hpp, this contains verbiage about H_DIV so we assign it the
1025  // FMS type FMS_HDIV.
1026 
1027  // I've seen RT_3D_P1 return the wrong order so get it from the name.
1028  int idim, iorder;
1029  if (sscanf(fecoll->Name(), "RT_%dD_P%d", &idim, &iorder) == 2)
1030  {
1031  order = (FmsInt)iorder;
1032  }
1033  else
1034  {
1035  order = static_cast<FmsInt>(fespace->GetOrder(0));
1036  }
1037 
1038  // Get the vdim from the fespace since the grid function is returning
1039  // 3 but we need it to be what was read from the file so we can pass the
1040  // right vdim to the FMS field descriptor to compute the expected number of dofs.
1041  vdim = fespace->GetVDim();
1042  break;
1043  }
1044  default:
1045  mfem::out << "Warning, unsupported ContType for field " << field_name
1046  << ". Using FMS_CONTINUOUS." << std::endl;
1047  ftype = FMS_CONTINUOUS;
1048  break;
1049  }
1050 
1051  // Now that we're not failing, create the fd and field.
1052  FmsFieldDescriptor fd = NULL;
1053  FmsField f = NULL;
1054  FmsDataCollectionAddFieldDescriptor(dc, fd_name.c_str(), &fd);
1055  FmsDataCollectionAddField(dc, field_name.c_str(), &f);
1056  *outfield = f;
1057 
1058  /* Q: Why is order defined on a per element basis? */
1059  FmsFieldDescriptorSetComponent(fd, comp);
1060  FmsFieldDescriptorSetFixedOrder(fd, ftype, btype, order);
1061 
1062  FmsInt ndofs;
1063  FmsFieldDescriptorGetNumDofs(fd, &ndofs);
1064 
1065  const char *name = NULL;
1066  FmsFieldGetName(f, &name);
1067  FmsLayoutType layout = fespace->GetOrdering() == mfem::Ordering::byVDIM ?
1068  FMS_BY_VDIM : FMS_BY_NODES;
1069 
1070 #ifdef DEBUG_MFEM_FMS
1071  switch (ftype)
1072  {
1073  case FMS_CONTINUOUS:
1074  mfem::out << "\tFMS_CONTINUOUS" << std::endl;
1075  break;
1076  case FMS_DISCONTINUOUS:
1077  mfem::out << "\tFMS_DISCONTINUOUS" << std::endl;
1078  break;
1079  case FMS_HDIV:
1080  mfem::out << "\tFMS_HDIV" << std::endl;
1081  break;
1082  }
1083  mfem::out << "\tField is order " << order << " with vdim " << vdim <<
1084  " and nDoFs " << ndofs << std::endl;
1085  mfem::out << "\tgf->size() " << gf->Size() << " ndofs * vdim " << ndofs * vdim
1086  << std::endl;
1087  mfem::out << "\tlayout " << layout << " (0 = BY_NODES, 1 = BY_VDIM)" <<
1088  std::endl;
1089 #endif
1090 
1091  if (FmsFieldSet(f, fd, vdim, layout, FMS_DOUBLE, c))
1092  {
1093  mfem::err << "Error setting field " << field_name << " in FMS." << std::endl;
1094  return 8;
1095  }
1096  return 0;
1097 }
1098 
1099 bool
1100 MfemMetaDataToFmsMetaData(DataCollection *mdc, FmsDataCollection fdc)
1101 {
1102  if (!mdc) { return false; }
1103  if (!fdc) { return false; }
1104 
1105  int *cycle = NULL;
1106  double *time = NULL, *timestep = NULL;
1107  FmsMetaData top_level = NULL;
1108  FmsMetaData *cycle_time_timestep = NULL;
1109  int mdata_err = 0;
1110  mdata_err = FmsDataCollectionAttachMetaData(fdc, &top_level);
1111  if (!top_level || mdata_err)
1112  {
1113  mfem::err << "Failed to attach metadata to the FmsDataCollection" << std::endl;
1114  return 40;
1115  }
1116 
1117  mdata_err = FmsMetaDataSetMetaData(top_level, "MetaData", 3,
1118  &cycle_time_timestep);
1119  if (!cycle_time_timestep || mdata_err)
1120  {
1121  mfem::err << "Failed to acquire FmsMetaData array" << std::endl;
1122  return false;
1123  }
1124 
1125  if (!cycle_time_timestep[0])
1126  {
1127  mfem::err << "The MetaData pointer for cycle is NULL" << std::endl;
1128  return false;
1129  }
1130  mdata_err = FmsMetaDataSetIntegers(cycle_time_timestep[0], "cycle",
1131  FMS_INT32, 1, (void**)&cycle);
1132  if (!cycle || mdata_err)
1133  {
1134  mfem::err << "The data pointer for cycle is NULL" << std::endl;
1135  return false;
1136  }
1137  *cycle = mdc->GetCycle();
1138 
1139  if (!cycle_time_timestep[1])
1140  {
1141  mfem::err << "The FmsMetaData pointer for time is NULL" << std::endl;
1142  return false;
1143  }
1144  mdata_err = FmsMetaDataSetScalars(cycle_time_timestep[1], "time", FMS_DOUBLE,
1145  1, (void**)&time);
1146  if (!time || mdata_err)
1147  {
1148  mfem::err << "The data pointer for time is NULL." << std::endl;
1149  return false;
1150  }
1151  *time = mdc->GetTime();
1152 
1153  if (!cycle_time_timestep[2])
1154  {
1155  mfem::err << "The FmsMetData pointer for timestep is NULL" << std::endl;
1156  return false;
1157  }
1158  mdata_err = FmsMetaDataSetScalars(cycle_time_timestep[2], "timestep",
1159  FMS_DOUBLE, 1, (void**)&timestep);
1160  if (!timestep || mdata_err)
1161  {
1162  mfem::err << "The data pointer for timestep is NULL" << std::endl;
1163  return false;
1164  }
1165  *timestep = mdc->GetTimeStep();
1166 
1167  return true;
1168 }
1169 
1170 //---------------------------------------------------------------------------
1171 bool
1172 FmsMetaDataGetInteger(FmsMetaData mdata, const std::string &key,
1173  std::vector<int> &values)
1174 {
1175  if (!mdata) { return false; }
1176 
1177  bool retval = false;
1178  FmsMetaDataType type;
1179  FmsIntType int_type;
1180  FmsInt i, size;
1181  FmsMetaData *children = nullptr;
1182  const void *data = nullptr;
1183  const char *mdata_name = nullptr;
1184  if (FmsMetaDataGetType(mdata, &type) == 0)
1185  {
1186  switch (type)
1187  {
1188  case FMS_INTEGER:
1189  if (FmsMetaDataGetIntegers(mdata, &mdata_name, &int_type, &size, &data) == 0)
1190  {
1191  if (strcasecmp(key.c_str(), mdata_name) == 0)
1192  {
1193  retval = true;
1194 
1195  // Interpret the integers and store them in the std::vector<int>
1196  switch (int_type)
1197  {
1198  case FMS_INT8:
1199  for (i = 0; i < size; i++)
1200  {
1201  values.push_back(static_cast<int>(reinterpret_cast<const int8_t*>(data)[i]));
1202  }
1203  break;
1204  case FMS_INT16:
1205  for (i = 0; i < size; i++)
1206  {
1207  values.push_back(static_cast<int>(reinterpret_cast<const int16_t*>(data)[i]));
1208  }
1209  break;
1210  case FMS_INT32:
1211  for (i = 0; i < size; i++)
1212  {
1213  values.push_back(static_cast<int>(reinterpret_cast<const int32_t*>(data)[i]));
1214  }
1215  break;
1216  case FMS_INT64:
1217  for (i = 0; i < size; i++)
1218  {
1219  values.push_back(static_cast<int>(reinterpret_cast<const int64_t*>(data)[i]));
1220  }
1221  break;
1222  case FMS_UINT8:
1223  for (i = 0; i < size; i++)
1224  {
1225  values.push_back(static_cast<int>(reinterpret_cast<const uint8_t*>(data)[i]));
1226  }
1227  break;
1228  case FMS_UINT16:
1229  for (i = 0; i < size; i++)
1230  {
1231  values.push_back(static_cast<int>(reinterpret_cast<const uint16_t*>(data)[i]));
1232  }
1233  break;
1234  case FMS_UINT32:
1235  for (i = 0; i < size; i++)
1236  {
1237  values.push_back(static_cast<int>(reinterpret_cast<const uint32_t*>(data)[i]));
1238  }
1239  break;
1240  case FMS_UINT64:
1241  for (i = 0; i < size; i++)
1242  {
1243  values.push_back(static_cast<int>(reinterpret_cast<const uint64_t*>(data)[i]));
1244  }
1245  break;
1246  default:
1247  retval = false;
1248  break;
1249  }
1250  }
1251  }
1252  break;
1253  case FMS_META_DATA:
1254  if (FmsMetaDataGetMetaData(mdata, &mdata_name, &size, &children) == 0)
1255  {
1256  // Recurse to look for the key we want.
1257  for (i = 0; i < size && !retval; i++)
1258  {
1259  retval = FmsMetaDataGetInteger(children[i], key, values);
1260  }
1261  }
1262  break;
1263  default:
1264  break;
1265  }
1266  }
1267 
1268  return retval;
1269 }
1270 
1271 //---------------------------------------------------------------------------
1272 bool
1273 FmsMetaDataGetScalar(FmsMetaData mdata, const std::string &key,
1274  std::vector<double> &values)
1275 {
1276  if (!mdata) { return false; }
1277 
1278  bool retval = false;
1279  FmsMetaDataType type;
1280  FmsScalarType scal_type;
1281  FmsInt i, size;
1282  FmsMetaData *children = nullptr;
1283  const void *data = nullptr;
1284  const char *mdata_name = nullptr;
1285  if (FmsMetaDataGetType(mdata, &type) == 0)
1286  {
1287  switch (type)
1288  {
1289  case FMS_SCALAR:
1290  if (FmsMetaDataGetScalars(mdata, &mdata_name, &scal_type, &size, &data) == 0)
1291  {
1292  if (strcasecmp(key.c_str(), mdata_name) == 0)
1293  {
1294  retval = true;
1295 
1296  // Interpret the integers and store them in the std::vector<int>
1297  switch (scal_type)
1298  {
1299  case FMS_FLOAT:
1300  for (i = 0; i < size; i++)
1301  {
1302  values.push_back(static_cast<double>(reinterpret_cast<const float*>(data)[i]));
1303  }
1304  break;
1305  case FMS_DOUBLE:
1306  for (i = 0; i < size; i++)
1307  {
1308  values.push_back(reinterpret_cast<const double*>(data)[i]);
1309  }
1310  break;
1311  default:
1312  retval = false;
1313  break;
1314  }
1315  }
1316  }
1317  break;
1318  case FMS_META_DATA:
1319  if (FmsMetaDataGetMetaData(mdata, &mdata_name, &size, &children) == 0)
1320  {
1321  // Recurse to look for the key we want.
1322  for (i = 0; i < size && !retval; i++)
1323  {
1324  retval = FmsMetaDataGetScalar(children[i], key, values);
1325  }
1326  }
1327  break;
1328  default:
1329  break;
1330  }
1331  }
1332 
1333  return retval;
1334 }
1335 
1336 //---------------------------------------------------------------------------
1337 bool
1338 FmsMetaDataGetString(FmsMetaData mdata, const std::string &key,
1339  std::string &value)
1340 {
1341  if (!mdata) { return false; }
1342 
1343  bool retval = false;
1344  FmsMetaDataType type;
1345  FmsInt i, size;
1346  FmsMetaData *children = nullptr;
1347  const char *mdata_name = nullptr;
1348  const char *str_value = nullptr;
1349 
1350  if (FmsMetaDataGetType(mdata, &type) == 0)
1351  {
1352  switch (type)
1353  {
1354  case FMS_STRING:
1355  if (FmsMetaDataGetString(mdata, &mdata_name, &str_value) == 0)
1356  {
1357  if (strcasecmp(key.c_str(), mdata_name) == 0)
1358  {
1359  retval = true;
1360  value = str_value;
1361  }
1362  }
1363  break;
1364  case FMS_META_DATA:
1365  if (FmsMetaDataGetMetaData(mdata, &mdata_name, &size, &children) == 0)
1366  {
1367  // Recurse to look for the key we want.
1368  for (i = 0; i < size && !retval; i++)
1369  {
1370  retval = FmsMetaDataGetString(children[i], key, value);
1371  }
1372  }
1373  break;
1374  default:
1375  break;
1376  }
1377  }
1378 
1379  return retval;
1380 }
1381 
1382 /* -------------------------------------------------------------------------- */
1383 /* FMS to MFEM conversion function */
1384 /* -------------------------------------------------------------------------- */
1385 
1386 int FmsDataCollectionToDataCollection(FmsDataCollection dc,
1387  DataCollection **mfem_dc)
1388 {
1389  int retval = 0;
1390  FmsMesh fms_mesh;
1391  FmsDataCollectionGetMesh(dc, &fms_mesh);
1392 
1393  // NOTE: The MFEM data collection has a single Mesh. Mesh has a constructor
1394  // to take multiple Mesh objects but it appears to glue them together.
1395  Mesh *mesh = nullptr;
1396  int err = FmsMeshToMesh(fms_mesh, &mesh);
1397  if (err == 0)
1398  {
1399  std::string collection_name("collection");
1400  const char *cn = nullptr;
1401  FmsDataCollectionGetName(dc, &cn);
1402  if (cn != nullptr)
1403  {
1404  collection_name = cn;
1405  }
1406 
1407  // Make a data collection that contains the mesh.
1408  DataCollection *mdc = new DataCollection(collection_name, mesh);
1409  mdc->SetOwnData(true);
1410 
1411  // Now do fields, etc. and add them to mdc.
1412  FmsField *fields = nullptr;
1413  FmsInt num_fields = 0;
1414  if (FmsDataCollectionGetFields(dc, &fields, &num_fields) == 0)
1415  {
1416  for (FmsInt i = 0; i < num_fields; ++i)
1417  {
1418  const char *name = nullptr;
1419  FmsFieldGetName(fields[i], &name);
1420 
1421  GridFunction *gf = new GridFunction;
1422 
1423  // Get the data type.
1424  FmsFieldDescriptor f_fd;
1425  FmsLayoutType f_layout;
1426  FmsScalarType f_data_type;
1427  const void *f_data;
1428  FmsFieldGet(fields[i], &f_fd, NULL, &f_layout, &f_data_type,
1429  &f_data);
1430 
1431  // Interpret the field according to its data type.
1432  int err = 1;
1433  switch (f_data_type)
1434  {
1435  case FMS_FLOAT:
1436  err = FmsFieldToGridFunction<float>(fms_mesh, fields[i], mesh, *gf, true);
1437  break;
1438  case FMS_DOUBLE:
1439  err = FmsFieldToGridFunction<double>(fms_mesh, fields[i], mesh, *gf, true);
1440  break;
1441  case FMS_COMPLEX_FLOAT:
1442  case FMS_COMPLEX_DOUBLE:
1443  // Does MFEM support complex?
1444  break;
1445  default:
1446  break;
1447  }
1448 
1449  if (err == 0)
1450  {
1451  mdc->RegisterField(name, gf);
1452  }
1453  else
1454  {
1455  mfem::out << "There was an error converting " << name << " code: " << err <<
1456  std::endl;
1457  delete gf;
1458  }
1459 
1460  const char *fname = NULL;
1461  FmsFieldGetName(fields[i], &fname);
1462  }
1463  }
1464 
1465  // If we have metadata in FMS, pass what we can through to MFEM.
1466  FmsMetaData mdata = NULL;
1467  FmsDataCollectionGetMetaData(dc, &mdata);
1468  if (mdata)
1469  {
1470  std::vector<int> ivalues;
1471  std::vector<double> dvalues;
1472  std::string svalue;
1473  if (FmsMetaDataGetInteger(mdata, "cycle", ivalues))
1474  {
1475  if (!ivalues.empty())
1476  {
1477  mdc->SetCycle(ivalues[0]);
1478  }
1479  }
1480  if (FmsMetaDataGetScalar(mdata, "time", dvalues))
1481  {
1482  if (!dvalues.empty())
1483  {
1484  mdc->SetTime(dvalues[0]);
1485  }
1486  }
1487  if (FmsMetaDataGetScalar(mdata, "timestep", dvalues))
1488  {
1489  if (!dvalues.empty())
1490  {
1491  mdc->SetTimeStep(dvalues[0]);
1492  }
1493  }
1494  }
1495 
1496  *mfem_dc = mdc;
1497  }
1498  else
1499  {
1500  mfem::out << "FmsDataCollectionToDataCollection: mesh failed to convert. err="
1501  << err
1502  << endl;
1503 
1504  retval = 1;
1505  }
1506 
1507 #if DEBUG_FMS_MFEM
1508  if (*mfem_dc)
1509  {
1510  VisItDataCollection visit_dc(std::string("DEBUG_DC"), mesh);
1511  visit_dc.SetOwnData(false);
1512  const auto &fields = (*mfem_dc)->GetFieldMap();
1513  for (const auto &field : fields)
1514  {
1515  visit_dc.RegisterField(field.first, field.second);
1516  }
1517  visit_dc.Save();
1518  }
1519 #endif
1520 
1521  return retval;
1522 }
1523 
1524 
1525 
1526 /* -------------------------------------------------------------------------- */
1527 /* MFEM to FMS conversion function */
1528 /* -------------------------------------------------------------------------- */
1529 
1530 int
1531 MeshToFmsMesh(const Mesh *mmesh, FmsMesh *fmesh, FmsComponent *volume)
1532 {
1533  if (!mmesh) { return 1; }
1534  if (!fmesh) { return 2; }
1535  if (!volume) { return 3; }
1536 
1537 
1538  int err = 0;
1539  const int num_verticies = mmesh->GetNV();
1540 #ifdef DEBUG_MFEM_FMS
1541  const int num_edges = mmesh->GetNEdges();
1542 #endif
1543  const int num_faces = mmesh->GetNFaces();
1544  const int num_elements = mmesh->GetNE();
1545 
1546 #ifdef DEBUG_MFEM_FMS
1547  mfem::out << "nverts: " << num_verticies << std::endl;
1548  mfem::out << "nedges: " << num_edges << std::endl;
1549  mfem::out << "nfaces: " << num_faces << std::endl;
1550  mfem::out << "nele: " << num_elements << std::endl;
1551 #endif
1552 
1553  FmsMeshConstruct(fmesh);
1554  FmsMeshSetPartitionId(*fmesh, 0, 1);
1555 
1556  FmsDomain *domains = NULL;
1557  FmsMeshAddDomains(*fmesh, "Domain", 1, &domains);
1558  FmsDomainSetNumVertices(domains[0], num_verticies);
1559 
1560  const int edge_reorder[2] = {1, 0};
1561  const int quad_reorder[4] = {0,1,2,3};
1562  const int tet_reorder[4] = {3,2,1,0};
1563  const int hex_reorder[6] = {0,5,1,3,4,2};
1564  const int *reorder[8] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL};
1565  reorder[FMS_EDGE] = edge_reorder;
1566  reorder[FMS_QUADRILATERAL] = quad_reorder;
1567  reorder[FMS_TETRAHEDRON] = tet_reorder;
1568  reorder[FMS_HEXAHEDRON] = hex_reorder;
1569 
1570  const mfem::Table *edges = mmesh->GetEdgeVertexTable();
1571  if (!edges)
1572  {
1573  mfem::err << "Error, mesh has no edges." << std::endl;
1574  return 1;
1575  }
1576  mfem::Table *faces = mmesh->GetFaceEdgeTable();
1577  if (!faces && num_faces > 0)
1578  {
1579  mfem::err <<
1580  "Error, mesh contains faces but the \"GetFaceEdgeTable\" returned NULL." <<
1581  std::endl;
1582  return 2;
1583  }
1584 
1585  // Build edges
1586  std::vector<int> edge_verts(edges->Size() * 2);
1587  for (int i = 0; i < edges->Size(); i++)
1588  {
1589  mfem::Array<int> nids;
1590  edges->GetRow(i, nids);
1591  for (int j = 0; j < 2; j++)
1592  {
1593  edge_verts[i*2 + j] = nids[j];
1594  }
1595  }
1596 
1597  // TODO: Move this code to after the for loop so edges can be added at top level entities
1598  FmsDomainSetNumEntities(domains[0], FMS_EDGE, FMS_INT32, edge_verts.size() / 2);
1599  FmsDomainAddEntities(domains[0], FMS_EDGE, reorder, FMS_INT32,
1600  edge_verts.data(), edge_verts.size() / 2);
1601 #ifdef DEBUG_MFEM_FMS
1602  mfem::out << "EDGES: ";
1603  for (int i = 0; i < edge_verts.size(); i++)
1604  {
1605  if (i % 2 == 0) { mfem::out << std::endl << "\t" << i/2 << ": "; }
1606  mfem::out << edge_verts[i] << " ";
1607  }
1608  mfem::out << std::endl;
1609 #endif
1610 
1611  // Build faces
1612  if (faces)
1613  {
1614  // TODO: Support Triangles and Quads, and move this code after the for
1615  // loop so these can be added as top level entities
1616  int rowsize = faces->RowSize(0);
1617  std::vector<int> face_edges(faces->Size() * rowsize);
1618  for (int i = 0; i < faces->Size(); i++)
1619  {
1620  mfem::Array<int> eids;
1621  faces->GetRow(i, eids);
1622  for (int j = 0; j < rowsize; j++)
1623  {
1624  face_edges[i*rowsize + j] = eids[j];
1625  }
1626  }
1627  FmsEntityType ent_type = (rowsize == 3) ? FMS_TRIANGLE : FMS_QUADRILATERAL;
1628  FmsDomainSetNumEntities(domains[0], ent_type, FMS_INT32,
1629  face_edges.size() / rowsize);
1630  FmsDomainAddEntities(domains[0], ent_type, NULL, FMS_INT32, face_edges.data(),
1631  face_edges.size() / rowsize);
1632 #ifdef DEBUG_MFEM_FMS
1633  mfem::out << "FACES: ";
1634  for (int i = 0; i < face_edges.size(); i++)
1635  {
1636  if (i % rowsize == 0) { mfem::out << std::endl << "\t" << i/rowsize << ": "; }
1637  mfem::out << "(" << edge_verts[face_edges[i]*2] << ", " <<
1638  edge_verts[face_edges[i]*2+1] << ") ";
1639  }
1640  mfem::out << std::endl;
1641 #endif
1642  }
1643 
1644  // Add top level elements
1645  std::vector<int> tags;
1646  std::vector<int> tris;
1647  std::vector<int> quads;
1648  std::vector<int> tets;
1649  std::vector<int> hexes;
1650  for (int i = 0; i < num_elements; i++)
1651  {
1652  auto etype = mmesh->GetElementType(i);
1653  tags.push_back(mmesh->GetAttribute(i));
1654  switch (etype)
1655  {
1656  case mfem::Element::POINT:
1657  {
1658  // TODO: ?
1659  break;
1660  }
1662  {
1663  // TODO: ?
1664  break;
1665  }
1667  {
1668  mfem::Array<int> eids, oris;
1669  mmesh->GetElementEdges(i, eids, oris);
1670  for (int e = 0; e < 3; e++)
1671  {
1672  tris.push_back(eids[e]);
1673  }
1674  break;
1675  }
1677  {
1678  mfem::Array<int> eids, oris;
1679  mmesh->GetElementEdges(i, eids, oris);
1680  for (int e = 0; e < 4; e++)
1681  {
1682  quads.push_back(eids[e]);
1683  }
1684  break;
1685  }
1687  {
1688  mfem::Array<int> fids, oris;
1689  mmesh->GetElementFaces(i, fids, oris);
1690  for (int f = 0; f < 4; f++)
1691  {
1692  tets.push_back(fids[f]);
1693  }
1694  break;
1695  }
1697  {
1698  mfem::Array<int> fids, oris;
1699  mmesh->GetElementFaces(i, fids, oris);
1700  for (int f = 0; f < 6; f++)
1701  {
1702  hexes.push_back(fids[f]);
1703  }
1704  break;
1705  }
1706  default:
1707  mfem::err << "Error, element not implemented." << std::endl;
1708  return 3;
1709  }
1710  }
1711 
1712  if (tris.size())
1713  {
1714  FmsDomainSetNumEntities(domains[0], FMS_TRIANGLE, FMS_INT32, tris.size() / 3);
1715  FmsDomainAddEntities(domains[0], FMS_TRIANGLE, reorder, FMS_INT32, tris.data(),
1716  tris.size() / 3);
1717 #ifdef DEBUG_MFEM_FMS
1718  mfem::out << "TRIS: ";
1719  for (int i = 0; i < tris.size(); i++)
1720  {
1721  if (i % 3 == 0) { mfem::out << std::endl << "\t" << i/3 << ": "; }
1722  mfem::out << tris[i] << " ";
1723  }
1724  mfem::out << std::endl;
1725 #endif
1726  }
1727 
1728  if (quads.size())
1729  {
1730  // TODO: Not quite right either, if there are hexes and quads then this
1731  // will overwrite the faces that made up the hexes
1732  FmsDomainSetNumEntities(domains[0], FMS_QUADRILATERAL, FMS_INT32,
1733  quads.size() / 4);
1734  FmsDomainAddEntities(domains[0], FMS_QUADRILATERAL, reorder, FMS_INT32,
1735  quads.data(), quads.size() / 4);
1736 #ifdef DEBUG_MFEM_FMS
1737  mfem::out << "QUADS: ";
1738  for (int i = 0; i < quads.size(); i++)
1739  {
1740  if (i % 4 == 0) { mfem::out << std::endl << "\t" << i/4 << ": "; }
1741  mfem::out << quads[i] << " ";
1742  }
1743  mfem::out << std::endl;
1744 #endif
1745  }
1746 
1747  if (tets.size())
1748  {
1749  FmsDomainSetNumEntities(domains[0], FMS_TETRAHEDRON, FMS_INT32,
1750  tets.size() / 4);
1751  FmsDomainAddEntities(domains[0], FMS_TETRAHEDRON, reorder, FMS_INT32,
1752  tets.data(), tets.size() / 4);
1753 #ifdef DEBUG_MFEM_FMS
1754  mfem::out << "TETS: ";
1755  for (int i = 0; i < tets.size(); i++)
1756  {
1757  if (i % 4 == 0) { mfem::out << std::endl << "\t" << i/4 << ": "; }
1758  mfem::out << tets[i] << " ";
1759  }
1760  mfem::out << std::endl;
1761 #endif
1762  }
1763 
1764  if (hexes.size())
1765  {
1766  FmsDomainSetNumEntities(domains[0], FMS_HEXAHEDRON, FMS_INT32,
1767  hexes.size() / 6);
1768  FmsDomainAddEntities(domains[0], FMS_HEXAHEDRON, reorder, FMS_INT32,
1769  hexes.data(), hexes.size() / 6);
1770 #ifdef DEBUG_MFEM_FMS
1771  mfem::out << "HEXES: ";
1772  for (int i = 0; i < hexes.size(); i++)
1773  {
1774  if (i % 6 == 0) { mfem::out << std::endl << "\t" << i/6 << ": "; }
1775  mfem::out << hexes[i] << " ";
1776  }
1777  mfem::out << std::endl;
1778 #endif
1779  }
1780 
1781  err = FmsMeshFinalize(*fmesh);
1782  if (err)
1783  {
1784  mfem::err << "FmsMeshFinalize returned error code " << err << std::endl;
1785  return 4;
1786  }
1787 
1788  err = FmsMeshValidate(*fmesh);
1789  if (err)
1790  {
1791  mfem::err << "FmsMeshValidate returned error code " << err << std::endl;
1792  return 5;
1793  }
1794 
1795  FmsComponent v = NULL;
1796  FmsMeshAddComponent(*fmesh, "volume", &v);
1797  FmsComponentAddDomain(v, domains[0]);
1798 
1799  FmsTag tag;
1800  FmsMeshAddTag(*fmesh, "element_attribute", &tag);
1801  FmsTagSetComponent(tag, v);
1802  FmsTagSet(tag, FMS_INT32, FMS_INT32, tags.data(), tags.size());
1803 
1804  // Add boundary component
1805  std::vector<int> bdr_eles[FMS_NUM_ENTITY_TYPES];
1806  std::vector<int> bdr_attributes;
1807  const int NBE = mmesh->GetNBE();
1808  for (int i = 0; i < NBE; i++)
1809  {
1810  const Element::Type betype = mmesh->GetBdrElementType(i);
1811  bdr_attributes.push_back(mmesh->GetBdrAttribute(i));
1812  switch (betype)
1813  {
1814  case Element::POINT:
1815  bdr_eles[FMS_VERTEX].push_back(mmesh->GetBdrElementEdgeIndex(i));
1816  break;
1817  case Element::SEGMENT:
1818  bdr_eles[FMS_EDGE].push_back(mmesh->GetBdrElementEdgeIndex(i));
1819  break;
1820  case Element::TRIANGLE:
1821  bdr_eles[FMS_TRIANGLE].push_back(mmesh->GetBdrElementEdgeIndex(i));
1822  break;
1824  bdr_eles[FMS_QUADRILATERAL].push_back(mmesh->GetBdrElementEdgeIndex(i));
1825  break;
1826  case Element::TETRAHEDRON:
1827  bdr_eles[FMS_TETRAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
1828  break;
1829  case Element::HEXAHEDRON:
1830  bdr_eles[FMS_HEXAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
1831  break;
1832  default:
1833  MFEM_WARNING("Unsupported boundary element " << betype << " at boundary index "
1834  << i);
1835  break;
1836  }
1837  }
1838 
1839  if (NBE)
1840  {
1841  FmsComponent boundary = NULL;
1842  FmsMeshAddComponent(*fmesh, "boundary", &boundary);
1843  FmsInt part_id;
1844  FmsComponentAddPart(boundary, domains[0], &part_id);
1845  for (int i = FMS_NUM_ENTITY_TYPES - 1; i > 0; i--)
1846  {
1847  if (bdr_eles[i].size())
1848  {
1849  FmsComponentAddPartEntities(boundary, part_id, (FmsEntityType)i,
1850  FMS_INT32, FMS_INT32, FMS_INT32, NULL,
1851  bdr_eles[i].data(),
1852  NULL, bdr_eles[i].size());
1853  break;
1854  }
1855  }
1856  FmsComponentAddRelation(v, 1);
1857  FmsTag boundary_tag = NULL;
1858  FmsMeshAddTag(*fmesh, "boundary_attribute", &boundary_tag);
1859  FmsTagSetComponent(boundary_tag, boundary);
1860  FmsTagSet(boundary_tag, FMS_INT32, FMS_INT32, bdr_attributes.data(),
1861  bdr_attributes.size());
1862  }
1863  *volume = v;
1864  return 0;
1865 }
1866 
1867 int
1869  FmsDataCollection *dc)
1870 {
1871  // TODO: Write me
1872  int err = 0;
1873  const Mesh *mmesh = mfem_dc->GetMesh();
1874 
1875  FmsMesh fmesh = NULL;
1876  FmsComponent volume = NULL;
1877  err = MeshToFmsMesh(mmesh, &fmesh, &volume);
1878  if (!fmesh || !volume || err)
1879  {
1880  mfem::err << "Error converting mesh topology from MFEM to FMS" << std::endl;
1881  if (fmesh)
1882  {
1883  FmsMeshDestroy(&fmesh);
1884  }
1885  return 1;
1886  }
1887 
1888  err = FmsDataCollectionCreate(fmesh, mfem_dc->GetCollectionName().c_str(), dc);
1889  if (!*dc || err)
1890  {
1891  mfem::err << "There was an error creating the FMS data collection." <<
1892  std::endl;
1893  FmsMeshDestroy(&fmesh);
1894  if (*dc)
1895  {
1896  FmsDataCollectionDestroy(dc);
1897  }
1898  return 3;
1899  }
1900 
1901  // Add the coordinates field to the data collection
1902  const mfem::GridFunction *mcoords = mmesh->GetNodes();
1903  if (mcoords)
1904  {
1905  FmsField fcoords = NULL;
1906  err = GridFunctionToFmsField(*dc, volume, "CoordsDescriptor", "Coords", mmesh,
1907  mcoords, &fcoords);
1908  err |= FmsComponentSetCoordinates(volume, fcoords);
1909  }
1910  else
1911  {
1912  // Sometimes the nodes are stored as just a vector of vertex coordinates
1913  mfem::Vector mverts;
1914  mmesh->GetVertices(mverts);
1915  FmsFieldDescriptor fdcoords = NULL;
1916  FmsField fcoords = NULL;
1917  err = FmsDataCollectionAddFieldDescriptor(*dc, "CoordsDescriptor", &fdcoords);
1918  err |= FmsFieldDescriptorSetComponent(fdcoords, volume);
1919  err |= FmsFieldDescriptorSetFixedOrder(fdcoords, FMS_CONTINUOUS,
1920  FMS_NODAL_GAUSS_CLOSED, 1);
1921  err |= FmsDataCollectionAddField(*dc, "Coords", &fcoords);
1922  err |= FmsFieldSet(fcoords, fdcoords, mmesh->SpaceDimension(), FMS_BY_NODES,
1923  FMS_DOUBLE, mverts.HostRead());
1924  err |= FmsComponentSetCoordinates(volume, fcoords);
1925  }
1926 
1927  if (err)
1928  {
1929  mfem::err << "There was an error setting the mesh coordinates." << std::endl;
1930  FmsMeshDestroy(&fmesh);
1931  FmsDataCollectionDestroy(dc);
1932  return 4;
1933  }
1934 
1935  const auto &fields = mfem_dc->GetFieldMap();
1936  for (const auto &pair : fields)
1937  {
1938  std::string fd_name(pair.first + "Descriptor");
1939  FmsField field;
1940  err = GridFunctionToFmsField(*dc, volume, fd_name, pair.first.c_str(), mmesh,
1941  pair.second, &field);
1942  if (err)
1943  {
1944  mfem::err << "WARNING: There was an error adding the " << pair.first <<
1945  " field. Continuing..." << std::endl;
1946  }
1947  }
1948 
1949  // /* TODO:
1950  // const auto &qfields = mfem_dc->GetQFieldMap();
1951  // for(const auto &pair : qfields) {
1952  // FmsFieldDescriptor fd = NULL;
1953  // FmsField f = NULL;
1954  // std::string fd_name(pair.first + "Collection");
1955  // FmsDataCollectionAddFieldDescriptor(*dc, fd_name.c_str(), &fd);
1956  // FmsDataCollectionAddField(*dc, pair.first.c_str(), &f);
1957  // GridFunctionToFmsField(*dc, fd, f, volume, pair.second); // TODO: Volume isn't always going to be correct
1958  // } */
1959 
1960  MfemMetaDataToFmsMetaData(mfem_dc, *dc);
1961 
1962  return 0;
1963 }
1964 
1965 } // namespace mfem
1966 #endif
bool FmsMetaDataGetInteger(FmsMetaData mdata, const std::string &key, std::vector< int > &values)
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
Definition: fespace.cpp:3130
Table * GetEdgeVertexTable() const
Definition: mesh.cpp:6419
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition: mesh.cpp:6298
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1694
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:6588
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
virtual int GetContType() const =0
bool FmsMetaDataGetScalar(FmsMetaData mdata, const std::string &key, std::vector< double > &values)
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
Normal component of vector field.
Definition: fe_coll.hpp:47
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
Definition: mesh.cpp:6649
bool MfemMetaDataToFmsMetaData(DataCollection *mdc, FmsDataCollection fdc)
T * GetData()
Returns the data.
Definition: array.hpp:115
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1680
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:1092
static const char * Name(int b_type)
Check and convert a BasisType constant to a string identifier.
Definition: fe_base.hpp:92
int GetCycle() const
Get time cycle (for time-dependent simulations)
int FmsMeshToMesh(FmsMesh fms_mesh, Mesh **mfem_mesh)
Definition: fmsconvert.cpp:452
int RowSize(int i) const
Definition: table.hpp:108
bool FmsMetaDataGetString(FmsMetaData mdata, const std::string &key, std::string &value)
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
virtual int DofForGeometry(Geometry::Type GeomType) const =0
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
int AddBdrTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1857
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1083
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1190
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1626
Nodes: x_i = (i+1/2)/n, i=0,...,n-1.
Definition: fe_base.hpp:36
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:6514
void GetVertices(Vector &vert_coord) const
Definition: mesh.cpp:8231
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
double b
Definition: lissajous.cpp:42
int AddTet(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1708
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:6382
Bernstein polynomials.
Definition: fe_base.hpp:33
double GetTimeStep() const
Get the simulation time step (for time-dependent simulations)
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition: mesh.hpp:1311
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
virtual const char * Name() const
Definition: fe_coll.hpp:80
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1843
int DataCollectionToFmsDataCollection(DataCollection *mfem_dc, FmsDataCollection *dc)
int FmsDataCollectionToDataCollection(FmsDataCollection dc, DataCollection **mfem_dc)
Closed GaussLegendre.
Definition: fe_base.hpp:38
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
void SetTime(double t)
Set physical time (for time-dependent simulations)
Tangential components of vector field.
Definition: fe_coll.hpp:46
void GetFaceInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified face.
Definition: fespace.cpp:3142
const FieldMapType & GetFieldMap() const
Get a const reference to the internal field map.
Table * GetFaceEdgeTable() const
Definition: mesh.cpp:6391
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
void SetOwnData(bool o)
Set the ownership of collection data.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1196
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
virtual void Save() override
Save the collection and a VisIt root file.
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Definition: mesh.cpp:1757
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2945
int GetOrder(int i) const
Returns the polynomial degree of the i&#39;th finite element.
Definition: fespace.hpp:692
int AddBdrQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1871
int GridFunctionToFmsField(FmsDataCollection dc, FmsComponent comp, const std::string &fd_name, const std::string &field_name, const Mesh *mesh, const GridFunction *gf, FmsField *outfield)
Definition: fmsconvert.cpp:948
Serendipity basis (squares / cubes)
Definition: fe_base.hpp:37
Nodes: x_i = (i+1)/(n+1), i=0,...,n-1.
Definition: fe_base.hpp:34
bool BasisTypeToFmsBasisType(int bt, FmsBasisType &btype)
Definition: fmsconvert.cpp:892
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
void SetTimeStep(double ts)
Set the simulation time step (for time-dependent simulations)
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:6644
int MeshToFmsMesh(const Mesh *mmesh, FmsMesh *fmesh, FmsComponent *volume)
int VectorDim() const
Definition: gridfunc.cpp:323
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:3042
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn&#39;t exist.
Definition: hash.hpp:700
int dim
Definition: ex24.cpp:53
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:195
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Mesh * GetMesh()
Get a pointer to the mesh in the collection.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:251
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe_base.hpp:35
Vector data type.
Definition: vector.hpp:58
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:1095
int GetId(int p1, int p2)
Get id of item whose parents are p1, p2... Create it if it doesn&#39;t exist.
Definition: hash.hpp:609
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
void GetElementInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified element.
Definition: fespace.cpp:3109
Field is continuous across element interfaces.
Definition: fe_coll.hpp:45
const std::string & GetCollectionName() const
Get the name of the collection.
int FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh, GridFunction &func, bool setFE)
This function converts an FmsField to an MFEM GridFunction.
Definition: fmsconvert.cpp:107
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
double GetTime() const
Get physical time (for time-dependent simulations)