Repository navigation
Incorrect visualization of submesh transferred to mesh #324
Description
Hello,
In the MFEM example below I have a simple 3D box geometry (geometry.zip) where I extract one of the faces as a 2D submesh. On this submesh, I create a vector grid function u_submesh and project some constant vector onto it. I then transfer u_submesh to the vector grid function u defined on the original 3D mesh.
When I use glvis to visualizeu_submesh on a 2D submesh, I see the correct values given by the projected vector.
However, when ask glvis to visualize u on the original 3D mesh, I see some unexpected values (see the figure at the end).
Meanwhile, by manually checking the u values on the corresponding face, I can verify that the transfer from u_submesh was correct.
Please help me understand why glvis shows wrong values in 3D case.
#include "mfem.hpp"
using namespace mfem;
using namespace std;
int main(int argc, char *argv[]) {
Mpi::Init(argc, argv);
Hypre::Init();
int order = 1; // Finite element polynomial degree
// 3D mesh -----------------------------------------------------------------
Mesh mesh_serial("netgen/geometry.nm");
ParMesh mesh(MPI_COMM_WORLD, mesh_serial);
mesh_serial.Clear(); // the serial mesh is no longer needed
mesh.Save("mesh"); // save mesh to the external file
ND_FECollection fec(order, mesh.Dimension());
ParFiniteElementSpace fespace(&mesh, &fec);
ParGridFunction u(&fespace);
u = 0;
// 2D submesh --------------------------------------------------------------
int iface = 4;
Array<int> marker(1);
marker[0] = iface + 1;
ParSubMesh submesh(ParSubMesh::CreateFromBoundary(mesh, marker));
submesh.Save("submesh");
ND_FECollection fec_submesh(order, submesh.Dimension());
ParFiniteElementSpace fespace_submesh(&submesh, &fec_submesh);
ParGridFunction u_submesh(&fespace_submesh);
u_submesh = 0;
// Initialize the submesh grid function ------------------------------------
Vector v(3);
v[0] = 1.0; v[1] = 2.0; v[2] = 3.0;
VectorConstantCoefficient vcoef(v);
u_submesh.ProjectCoefficient(vcoef);
u_submesh.Save("u_submesh");
// Transfer from submesh to mesh -------------------------------------------
submesh.Transfer(u_submesh, u);
u.Save("u");
// Inspect the boundary values of the grid function ------------------------
VectorGridFunctionCoefficient u_coef(&u);
Geometry geometry;
for(int bfi = 0; bfi < mesh.GetNBE(); ++bfi) {
const int bface_attr = mesh.GetBdrAttribute(bfi);
if(bface_attr != iface+1) {
continue;
}
const FiniteElement *fe = fespace.GetBE(bfi);
ElementTransformation *tr = mesh.GetBdrElementTransformation(bfi);
// Apply the transformation:
const auto ip = geometry.GetCenter(fe->GetGeomType());
tr->SetIntPoint(&ip);
Vector vec(3);
u_coef.Eval(vec, *tr, ip);
cout << "Boundary Vec: [" << vec(0) << ", "
<< vec(1) << ", "
<< vec(2) << "]" << endl;
}
return 0;
}The x component of u is expected to be 1 on the front face, but what I see is
Activity