Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incorrect visualization of submesh transferred to mesh #324

Open
vladimir-yu-fedorov opened this issue Dec 20, 2024 · 0 comments
Open

Incorrect visualization of submesh transferred to mesh #324

vladimir-yu-fedorov opened this issue Dec 20, 2024 · 0 comments
Assignees
Labels

Comments

@vladimir-yu-fedorov
Copy link

vladimir-yu-fedorov commented Dec 20, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

4 participants