Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <fstream>
- #include <cassert>
- #include "mfem.hpp"
- extern "C" {
- #include "xf.h"
- #include "xf_All.h"
- #include "xf_Mesh.h"
- #include "xf_Basis.h"
- #include "xf_Data.h"
- #include "xf_State.h"
- #include "xf_Math.h"
- }
- // Convert XFlow shapes to MFEM geometry
- mfem::Geometry::Type ShapeToGeometry(enum xfe_ShapeType shape)
- {
- switch (shape)
- {
- case xfe_Point : return mfem::Geometry::POINT;
- case xfe_Segment : return mfem::Geometry::SEGMENT;
- case xfe_Triangle : return mfem::Geometry::TRIANGLE;
- case xfe_Quadrilateral : return mfem::Geometry::SQUARE;
- case xfe_Hexahedron : return mfem::Geometry::CUBE;
- case xfe_Tetrahedron : return mfem::Geometry::TETRAHEDRON;
- default: throw std::runtime_error("Unknown shape");
- }
- }
- void ReverseIntArray(const int *input, int size, int *output)
- {
- for (int i = 0; i < size; i++) output[input[i]] = i;
- }
- mfem::GridFunction convertXFVector(const xf_Mesh *xflow_mesh,
- mfem::Mesh *mesh,
- const xf_Vector *V)
- {
- int ierr;
- const int vdim = V->StateRank;
- const int dim = xflow_mesh->Dim;
- // Make sure the vector is interpolated
- if (V->Order == NULL) throw std::runtime_error("Data is not interpolated");
- // ... and has element linkage
- if (V->Linkage != xfe_LinkageGlobElem) throw std::runtime_error("Is not implemented for non-element-linked data");
- // ... and that the number of element groups equals the number of arrays
- if (V->nArray != xflow_mesh->nElemGroup) throw std::runtime_error("Need number of elements groups to match number of arrays");
- // Find max order
- int max_order = 0;
- for (int array = 0; array < V->nArray; array++)
- {
- if (V->vOrder == NULL)
- {
- const int order = V->Order[array];
- max_order = order > max_order ? order : max_order;
- }
- else
- {
- for (int comp = 0; comp < V->nComp[array]; comp++)
- {
- const int order = V->vOrder[array][comp];
- max_order = order > max_order ? order : max_order;
- }
- }
- }
- enum xfe_ShapeType elem_shape;
- ierr = xf_Error(xf_Basis2Shape(xflow_mesh->ElemGroup[0].QBasis, &elem_shape));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- // elem_basis is the basis of the solution _not_ necessarily that of the mesh
- enum xfe_BasisType elem_basis;
- ierr = xf_Error(xf_Shape2UniformLagrange(elem_shape, &elem_basis));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- int nn;
- ierr = xf_Error(xf_Order2nNode(elem_basis, max_order, &nn));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- mfem::L2_FECollection *fec = new mfem::L2_FECollection(max_order, xflow_mesh->Dim, mfem::BasisType::ClosedUniform);
- mfem::FiniteElementSpace *fes = new mfem::FiniteElementSpace(mesh, fec, vdim, mfem::Ordering::byVDIM);
- mfem::GridFunction gf(fes);
- // TODO make sure fec gets deleted somehow...
- // gf.MakeOwner(fec); // gf will destroy 'fec' and 'fes'
- // Create a temporary data set
- xf_DataSet *DataSet;
- ierr = xf_Error(xf_CreateDataSet(&DataSet));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- // This is the same operation as filling the mesh Nodes
- // GridFunction, but it cannot be reused since xflow does not expose
- // these as an xf_Vector.
- mfem::Array<int> vdofs;
- xf_Matrix *T = NULL;
- for (int egrp = 0; egrp < V->nArray; egrp++)
- {
- xf_GenArray &GA = V->GenArray[egrp];
- for (int elem = 0; elem < GA.n; elem++)
- {
- const enum xfe_BasisType basis = (V->Basis != NULL) ? V->Basis[egrp] : V->vBasis[egrp][elem];
- const int order = (V->vOrder != NULL) ? V->vOrder[egrp][elem] : V->Order[egrp];
- fes->GetElementVDofs(elem, vdofs);
- const mfem::FiniteElement *fe = fes->GetFE(elem);
- const int ndof = fe->GetDof();
- // Assumption: element data (so, use elem_shape from above)
- mfem::Array<double> values;
- if (order != max_order)
- {
- // Variable order - need to find a transfer matrix in XFlow
- ierr = xf_Error(xf_FindTransferMatrix(DataSet, basis, order, elem_basis, max_order, &T));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- int ndof_e;
- ierr = xf_Error(xf_Order2nNode(basis, order, &ndof_e));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- values.SetSize(ndof_e * ndof);
- xf_MxM_Set(T->GenArray[0].rValue[0], GA.rValue[elem], ndof, ndof_e, vdim, values);
- }
- else
- {
- values.MakeRef(GA.rValue[elem], ndof * vdim);
- }
- // The ordering of MFEM grid functions is _almost_ the same as
- // XFlow (both store lexicographically) but MFEM stores vdofs
- // byNODES so we need to loop over the dimensions manually:
- for (int i = 0; i < ndof; i++) {
- for (int vd = 0; vd < vdim; vd++) {
- gf(vdofs[i+ndof*vd]) = values[i * vdim + vd];
- // This is the syntax if we had to map:
- // gf(dofs[dof_map[i]]*vdim + vd) = values[i * dim + vd];
- }
- }
- // If the ordering where the same we could simply use this:
- // gf.SetSubVector(vdofs, values);
- }
- }
- ierr = xf_Error(xf_DestroyDataSet(DataSet));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- return gf;
- }
- class NodeHash
- {
- private:
- std::vector<int> forward;
- std::vector<int> reverse; // indexed by node
- public:
- NodeHash() : forward(), reverse() { }
- int add(int node) // Returns the tag
- {
- if (reverse.size() < node+1)
- {
- reverse.resize(node+1, -1);
- }
- int tag = reverse[node];
- if (tag < 0)
- {
- reverse[node] = forward.size();
- forward.push_back(node);
- }
- return tag;
- }
- int get(int tag) const
- {
- return forward[tag];
- }
- int size() const
- {
- return forward.size();
- }
- };
- class XFMesh : public mfem::Mesh
- {
- public:
- XFMesh(xf_Mesh &xflow_mesh);
- };
- XFMesh::XFMesh(xf_Mesh &xflow_mesh) : mfem::Mesh()
- {
- int ierr;
- std::vector<int> dof_map_rev, verts, fvec;
- // Count elements
- int nElemTot = 0;
- for (int egrp = 0; egrp < xflow_mesh.nElemGroup; egrp++)
- {
- nElemTot += xflow_mesh.ElemGroup[egrp].nElem;
- }
- // Count boundary elements
- int nBdrElem = 0;
- for (int bfgrp = 0; bfgrp < xflow_mesh.nBFaceGroup; bfgrp++)
- {
- nBdrElem += xflow_mesh.BFaceGroup[bfgrp].nBFace;
- }
- const int dim = xflow_mesh.Dim;
- // Check whether or not QOrder changes. If so, we cannot
- // represent the mesh nodes by a simple H1 grid function.
- int pQOrder = xflow_mesh.ElemGroup[0].QOrder;
- bool QOrderChanges = false;
- int MaxQOrder = pQOrder;
- for (int egrp = 1; egrp < xflow_mesh.nElemGroup; egrp++)
- {
- QOrderChanges |= (xflow_mesh.ElemGroup[egrp].QOrder != pQOrder);
- MaxQOrder = max(xflow_mesh.ElemGroup[egrp].QOrder, MaxQOrder);
- }
- if (QOrderChanges)
- {
- std::cout << "QOrder changes - will not curve mesh. In the future this function could determine the high order nodes on the low order elements." << std::endl;
- }
- // Create a FiniteElementCollection early on for dof_map access
- // dof_map = Cartesian to local H1 dof map
- mfem::H1_FECollection *fec = new mfem::H1_FECollection(MaxQOrder, dim);
- // Node hash to get the vertex ordering
- NodeHash nh;
- for (int egrp = 0; egrp < xflow_mesh.nElemGroup; egrp++)
- {
- const xf_ElemGroup &EG = xflow_mesh.ElemGroup[egrp];
- enum xfe_ShapeType shape;
- ierr = xf_Basis2Shape(EG.QBasis, &shape);
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- const int geom_type = ShapeToGeometry(shape);
- const int *dof_map = fec->GetDofMap(geom_type);
- const mfem::FiniteElement *fe = fec->FiniteElementForGeometry(geom_type);
- const int num_dof = fe->GetDof();
- dof_map_rev.resize(num_dof);
- ReverseIntArray(dof_map, num_dof, dof_map_rev.data());
- const int num_verts = mfem::Geometry::NumVerts[geom_type];
- for (int elem = 0; elem < EG.nElem; elem++)
- {
- for (int i = 0; i < num_verts; i++)
- {
- nh.add(EG.Node[elem][dof_map_rev[i]]);
- }
- }
- }
- InitMesh(dim, dim, nh.size(), nElemTot, nBdrElem);
- // 1. Add the vertices
- for (int i = 0; i < nh.size(); i++)
- {
- const int node = nh.get(i);
- AddVertex(xflow_mesh.Coord[node]);
- }
- NumOfVertices = nh.size();
- // 2. Add the elements
- int egrp_offset = 0;
- for (int egrp = 0; egrp < xflow_mesh.nElemGroup; egrp++)
- {
- const xf_ElemGroup &EG = xflow_mesh.ElemGroup[egrp];
- mfem::Element **elements_egrp = elements.GetData() + egrp_offset;
- enum xfe_ShapeType shape;
- ierr = xf_Basis2Shape(EG.QBasis, &shape);
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- switch (shape)
- {
- case xfe_Segment:
- for (int e = 0; e < EG.nElem; e++) elements_egrp[e] = new mfem::Segment();
- break;
- case xfe_Quadrilateral:
- for (int e = 0; e < EG.nElem; e++) elements_egrp[e] = new mfem::Quadrilateral();
- break;
- case xfe_Triangle:
- for (int e = 0; e < EG.nElem; e++) elements_egrp[e] = new mfem::Triangle();
- break;
- case xfe_Hexahedron:
- for (int e = 0; e < EG.nElem; e++) elements_egrp[e] = new mfem::Hexahedron();
- break;
- case xfe_Tetrahedron:
- for (int e = 0; e < EG.nElem; e++) elements_egrp[e] = new mfem::Tetrahedron();
- break;
- }
- const int geom_type = ShapeToGeometry(shape);
- const int *dof_map = fec->GetDofMap(geom_type);
- const mfem::FiniteElement *fe = fec->FiniteElementForGeometry(geom_type);
- const int num_dof = fe->GetDof();
- dof_map_rev.resize(num_dof);
- ReverseIntArray(dof_map, num_dof, dof_map_rev.data());
- const int num_verts = mfem::Geometry::NumVerts[geom_type];
- verts.resize(num_verts);
- for (int elem = 0; elem < EG.nElem; elem++)
- {
- for (int i = 0; i < num_verts; i++)
- {
- verts[i] = nh.add(EG.Node[elem][dof_map_rev[i]]);
- }
- elements_egrp[elem]->SetVertices(verts.data());
- elements_egrp[elem]->SetAttribute(1 + egrp);
- }
- egrp_offset += EG.nElem;
- }
- NumOfElements = egrp_offset;
- // 3. Add the boundary elements
- int bfgrp_offset = 0;
- for (int bfgrp = 0; bfgrp < xflow_mesh.nBFaceGroup; bfgrp++)
- {
- xf_BFaceGroup &BFG = xflow_mesh.BFaceGroup[bfgrp];
- mfem::Element **boundary_bfgrp = boundary.GetData() + bfgrp_offset;
- for (int bface = 0; bface < BFG.nBFace; bface++)
- {
- // Every element in the boundary face group may have a different
- // shape or derive from an element of different QOrder, hence
- // the allocation needs to happen inside the loop over bface.
- const xf_BFace &BF = BFG.BFace[bface];
- const xf_ElemGroup &EG = xflow_mesh.ElemGroup[BF.ElemGroup];
- enum xfe_ShapeType elem_shape, face_shape;
- ierr = xf_Basis2Shape(EG.QBasis, &elem_shape);
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- ierr = xf_FaceShape(elem_shape, BF.Face, &face_shape);
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- switch (face_shape)
- {
- case xfe_Point:
- boundary_bfgrp[bface] = new mfem::Point();
- break;
- case xfe_Segment:
- boundary_bfgrp[bface] = new mfem::Segment();
- break;
- case xfe_Quadrilateral:
- boundary_bfgrp[bface] = new mfem::Quadrilateral();
- break;
- case xfe_Triangle:
- boundary_bfgrp[bface] = new mfem::Triangle();
- break;
- }
- const int geom_type = ShapeToGeometry(face_shape);
- const mfem::FiniteElement *fe = fec->FiniteElementForGeometry(geom_type);
- int num_dof;
- if (geom_type != mfem::Geometry::POINT)
- {
- const int *dof_map = fec->GetDofMap(geom_type);
- num_dof = fe->GetDof();
- ReverseIntArray(dof_map, num_dof, dof_map_rev.data());
- }
- else
- {
- const int dof_map[1] = {0};
- num_dof = 1;
- ReverseIntArray(dof_map, num_dof, dof_map_rev.data());
- }
- dof_map_rev.resize(num_dof);
- const int num_verts = mfem::Geometry::NumVerts[geom_type];
- int nqnode;
- fvec.resize(num_dof);
- ierr = xf_NodesOnFace(EG.QBasis, EG.QOrder, BF.Face, &nqnode, fvec.data());
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- verts.resize(num_verts);
- for (int i = 0; i < num_verts; i++)
- {
- verts[i] = nh.add(EG.Node[BF.Elem][fvec[dof_map_rev[i]]]);
- }
- boundary_bfgrp[bface]->SetVertices(verts.data());
- boundary_bfgrp[bface]->SetAttribute(1 + bfgrp);
- }
- bfgrp_offset += BFG.nBFace;
- }
- NumOfBdrElements = bfgrp_offset;
- FinalizeTopology();
- if (MaxQOrder == 1) return;
- // --- Set curvature to MaxQOrder ---
- // The method below works, but it does extra work
- // SetCurvature(MaxQOrder, false, dim, mfem::Ordering::byVDIM);
- // Adding back the Q1 nodes manually is much more efficient
- mfem::FiniteElementSpace *fes = new mfem::FiniteElementSpace(this, fec, dim,
- mfem::Ordering::byVDIM);
- Nodes = new mfem::GridFunction(fes);
- Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
- own_nodes = 1;
- mfem::Array<int> dofs;
- for (int egrp = 0; egrp < xflow_mesh.nElemGroup; egrp++)
- {
- xf_ElemGroup &EG = xflow_mesh.ElemGroup[egrp];
- enum xfe_ShapeType shape;
- ierr = xf_Basis2Shape(EG.QBasis, &shape);
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- const int geom_type = ShapeToGeometry(shape);
- const int *dof_map = fec->GetDofMap(geom_type);
- for (int elem = 0; elem < EG.nElem; elem++)
- {
- fes->GetElementDofs(elem, dofs);
- const mfem::FiniteElement *fe = fes->GetFE(elem);
- const int ndof = fe->GetDof();
- for (int i = 0; i < ndof; i++) {
- for (int vd = 0; vd < dim; vd++) {
- (*Nodes)(dofs[dof_map[i]]*dim + vd) = xflow_mesh.Coord[0][EG.Node[elem][i] * dim + vd];
- }
- }
- }
- }
- }
- int main(int argc, char *argv[])
- {
- int ierr;
- xf_All *All;
- const char *xf_mesh_file = "NULL";
- const char *xf_data_file = "NULL";
- const char *xf_xfa_file = "NULL";
- const char *convert_list_char = "";
- bool write = false;
- bool plot = true;
- /// GLVis options
- char vishost[] = "localhost";
- int visport = 19916;
- mfem::OptionsParser args(argc, argv);
- args.AddOption(&xf_mesh_file, "-m", "--mesh",
- "Mesh file to use.");
- args.AddOption(&xf_data_file, "-d", "--data",
- "Data file to use.");
- args.AddOption(&xf_xfa_file, "-x", "--xfa",
- "Xfa file to use.");
- args.AddOption(&convert_list_char, "-c", "--convert",
- "Comma-separated list of data names to convert.");
- args.AddOption(&write, "-w", "--write", "-no-w", "--no-write",
- "Enable or disable output writing.");
- args.AddOption(&plot, "-p", "--plot", "-no-p", "--no-plot",
- "Enable or disable plotting.");
- args.Parse();
- if (!args.Good())
- {
- args.PrintUsage(std::cout);
- return 1;
- }
- args.PrintOptions(std::cout);
- /* Create .xfa structure */
- xf_Call(xf_CreateAll(&All, xfe_False));
- if (strcmp(xf_mesh_file, "NULL"))
- {
- // Mesh file passed
- std::cout << "Using mesh " << xf_mesh_file << std::endl;
- xf_Call(xf_ReadGriFile(xf_mesh_file, NULL, All->Mesh));
- }
- else if (strcmp(xf_xfa_file, "NULL"))
- {
- // xfa file passed
- std::cout << "Using xfa " << xf_xfa_file << std::endl;
- xf_Call(xf_ReadAllBinary(xf_xfa_file, All));
- if (!strcmp(xf_data_file, "NULL"))
- {
- std::cout << "with " << xf_xfa_file;
- ierr = xf_Error(xf_ReadAllBinary(xf_xfa_file, All));
- if (ierr != xf_OK) throw std::runtime_error("An error occured");
- }
- }
- if (strcmp(xf_data_file, "NULL"))
- {
- // Additional data file passed
- std::cout << "Using data file " << xf_data_file << std::endl;
- xf_DataSet *DataSet;
- xf_Call(xf_CreateDataSet(&DataSet));
- xf_Call(xf_ReadDataSetBinary(All->Mesh, NULL, xf_data_file, DataSet));
- xf_Call(xf_DataSetMerge(DataSet, All->DataSet));
- xf_Call(xf_DestroyDataSet(DataSet));
- }
- std::cout << "Converting mesh" << std::flush;
- // --- Convert mesh ---
- XFMesh mesh(*All->Mesh);
- std::cout << " ... done" << std::endl;
- std::cout << "Reading data" << std::flush;
- std::string convert_list(convert_list_char);
- int num_converted = 0;
- std::vector<mfem::GridFunction> gf_list;
- std::vector<std::string> gf_names;
- for (xf_Data *D = All->DataSet->Head; D != NULL; D = D->Next)
- {
- std::string data_title(D->Title);
- const std::size_t pos = convert_list.find(data_title);
- if (pos == std::string::npos) continue;
- std::cout << "Converting " << data_title << std::endl;
- if (D->Type == xfe_Vector)
- {
- xf_Vector *V = (xf_Vector *) D->Data;
- if (V->Linkage == xfe_LinkageGlobElem)
- {
- gf_list.emplace_back(convertXFVector(All->Mesh, &mesh, V));
- gf_names.emplace_back(std::string(D->Title));
- num_converted++;
- }
- }
- else if (D->Type == xfe_VectorGroup)
- {
- xf_VectorGroup *VG = (xf_VectorGroup *) D->Data;
- for (int i = 0; i < VG->nVector; i++)
- {
- if (VG->Vector[i]->Linkage == xfe_LinkageGlobElem)
- {
- gf_list.emplace_back(convertXFVector(All->Mesh, &mesh, VG->Vector[i]));
- gf_names.emplace_back(std::string(D->Title) + "_Vector" + std::to_string(i));
- num_converted++;
- }
- }
- }
- else
- {
- mfem::mfem_error("Type not supported");
- }
- }
- std::cout << " ... done - number converted = " << num_converted << std::endl;
- if (write)
- {
- // Save mesh
- std::ofstream ofs("xflow.mesh");
- ofs.precision(8);
- mesh.Print(ofs);
- // Save gridfunctions
- for (int i = 0; i < gf_list.size(); i++)
- {
- std::stringstream ss;
- ss << gf_names[i] << ".gf";
- std::ofstream ofs(ss.str());
- ofs.precision(8);
- gf_list[i].Save(ofs);
- }
- }
- if (plot)
- {
- if (gf_list.size() == 0)
- {
- // Plot only the mesh
- mfem::socketstream sout(vishost, visport);
- sout.precision(8);
- sout << "mesh\n" << mesh << std::flush;
- }
- else
- {
- // Loop over gridfunctions
- int Wx = 0, Wy = 0; // window position
- const int Ww = 350, Wh = 350; // window size
- int offx = Ww+10; // window offsets
- for (int i = 0; i < gf_list.size(); i++)
- {
- mfem::socketstream sout(vishost, visport);
- if (gf_list[i].VectorDim() == 2)
- {
- sout << "solution\n" << mesh << gf_list[i] << std::flush;
- sout << "window_title '" << gf_names[i] << "'\n"
- << "window_geometry "
- << Wx << " " << Wy << " " << Ww << " " << Wh << "\n" << std::flush;
- }
- else
- {
- std::cout << "WARNING: Skipping " << gf_names[i] << " because it has more than one vector output. Use -write and plot manually with GLVis using the -gc option." << std::endl;
- }
- }
- }
- }
- /* Destroy .xfa structure */
- xf_Call(xf_DestroyAll(All));
- return 0;
- }
Add Comment
Please, Sign In to add comment