CheckXmlFile.cpp File Reference

#include <MultiRegions/ExpList.h>
#include <SpatialDomains/MeshGraphIO.h>

Go to the source code of this file.

Classes
struct	Ord

Functions
bool	CheckTetRotation (Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, Array< OneD, NekDouble > &xz, SpatialDomains::TetGeom *tet, std::map< int, int > &vid2cnt)
int	main (int argc, char *argv[])

Function Documentation

CheckTetRotation()

bool CheckTetRotation	(	Array< OneD, NekDouble > &	xc,
		Array< OneD, NekDouble > &	yc,
		Array< OneD, NekDouble > &	xz,
		SpatialDomains::TetGeom *	tet,
		std::map< int, int > &	vid2cnt
	)

Definition at line 216 of file CheckXmlFile.cpp.

{
    bool RotationOK = true;
    std::array<Ord, 4> v;
    NekDouble abx, aby, abz;
 
    v[0].x = xc[vid2cnt[tet->GetVid(0)]];
    v[0].y = yc[vid2cnt[tet->GetVid(0)]];
    v[0].z = zc[vid2cnt[tet->GetVid(0)]];
 
    v[1].x = xc[vid2cnt[tet->GetVid(1)]];
    v[1].y = yc[vid2cnt[tet->GetVid(1)]];
    v[1].z = zc[vid2cnt[tet->GetVid(1)]];
 
    v[2].x = xc[vid2cnt[tet->GetVid(2)]];
    v[2].y = yc[vid2cnt[tet->GetVid(2)]];
    v[2].z = zc[vid2cnt[tet->GetVid(2)]];
 
    v[3].x = xc[vid2cnt[tet->GetVid(3)]];
    v[3].y = yc[vid2cnt[tet->GetVid(3)]];
    v[3].z = zc[vid2cnt[tet->GetVid(3)]];
 
    // cross product of edge 0 and 2
    abx = (v[1].y - v[0].y) * (v[2].z - v[0].z) -
          (v[1].z - v[0].z) * (v[2].y - v[0].y);
    aby = (v[1].z - v[0].z) * (v[2].x - v[0].x) -
          (v[1].x - v[0].x) * (v[2].z - v[0].z);
    abz = (v[1].x - v[0].x) * (v[2].y - v[0].y) -
          (v[1].y - v[0].y) * (v[2].x - v[0].x);
 
    // inner product of cross product with respect to edge 3 should be positive
    if (((v[3].x - v[0].x) * abx + (v[3].y - v[0].y) * aby +
         (v[3].z - v[0].z) * abz) < 0.0)
    {
        std::cerr << "ERROR: Element " << tet->GetGlobalID()
                  << "is NOT counter-clockwise\n"
                  << std::endl;
        RotationOK = false;
    }
    return RotationOK;
}

References Nektar::SpatialDomains::Geometry::GetGlobalID(), and Nektar::SpatialDomains::Geometry::GetVid().

Referenced by main().

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 44 of file CheckXmlFile.cpp.

{
    Array<OneD, NekDouble> fce;
    Array<OneD, NekDouble> xc0, xc1, xc2;
 
    if (argc != 2)
    {
        std::cerr << "Usage: CheckXmlFile  meshfile.xml" << std::endl;
        return 1;
    }
 
    LibUtilities::SessionReaderSharedPtr vSession =
        LibUtilities::SessionReader::CreateInstance(argc, argv);
 
    //----------------------------------------------
    // Read in mesh from input file
    std::string meshfile(argv[argc - 1]);
    SpatialDomains::MeshGraphSharedPtr mesh =
        SpatialDomains::MeshGraphIO::Read(vSession);
    //----------------------------------------------
 
    //----------------------------------------------
    // Define Expansion
    int expdim = mesh->GetMeshDimension();
 
    switch (expdim)
    {
        case 1:
            ASSERTL0(false, "1D not set up");
            break;
        case 2:
            ASSERTL0(false, "2D not set up");
            break;
        case 3:
        {
            int cnt = 0, nverts = mesh->GetNvertices();
            Array<OneD, NekDouble> xc(nverts), yc(nverts), zc(nverts);
            std::map<int, int> vid2cnt;
 
            for (auto [id, vert] :
                 mesh->GetGeomMap<SpatialDomains::PointGeom>())
            {
                // Check element rotation
                vert->GetCoords(xc[cnt], yc[cnt], zc[cnt]);
                vid2cnt[vert->GetGlobalID()] = cnt++;
            }
 
            // Check for any duplicated vertex coordinates.
            for (int i = 0; i < nverts; ++i)
            {
                for (int j = i + 1; j < nverts; ++j)
                {
                    if ((xc[i] - xc[j]) * (xc[i] - xc[j]) +
                            (yc[i] - yc[j]) * (yc[i] - yc[j]) +
                            (zc[i] - zc[j]) * (zc[i] - zc[j]) <
                        1e-10)
                    {
                        std::cerr << "ERROR: Duplicate vertices: " << i << " "
                                  << j << std::endl;
                    }
                }
            }
 
            bool NoRotateIssues      = true;
            bool NoOrientationIssues = true;
            for (auto [id, tet] : mesh->GetGeomMap<SpatialDomains::TetGeom>())
            {
                // check rotation and dump
                NoOrientationIssues =
                    CheckTetRotation(xc, yc, zc, tet, vid2cnt);
 
                // Check face rotation
                if (tet->GetFace(0)->GetVid(2) != tet->GetVid(2))
                {
                    std::cerr << "ERROR: Face " << tet->GetFid(0) << " (vert "
                              << tet->GetFace(0)->GetVid(2)
                              << ") is not aligned with base vertex of Tet "
                              << tet->GetGlobalID() << " (vert "
                              << tet->GetVid(2) << ")" << std::endl;
                    NoRotateIssues = false;
                }
 
                for (int i = 1; i < 4; ++i)
                {
                    if (tet->GetFace(i)->GetVid(2) != tet->GetVid(3))
                    {
                        std::cerr << "ERROR: Face " << tet->GetFid(i)
                                  << " is not aligned with top Vertex of Tet "
                                  << tet->GetGlobalID() << std::endl;
                        NoRotateIssues = false;
                    }
                }
            }
            if (NoOrientationIssues)
            {
                std::cout << "All Tet have correct ordering for anticlockwise "
                             "rotation"
                          << std::endl;
            }
 
            if (NoRotateIssues)
            {
                std::cout << "All Tet faces are correctly aligned" << std::endl;
            }
 
            for ([[maybe_unused]] auto [id, pyr] :
                 mesh->GetGeomMap<SpatialDomains::PyrGeom>())
            {
                // Put pyramid checks in here
            }
 
            NoRotateIssues      = true;
            NoOrientationIssues = true;
            for (auto [id, pri] : mesh->GetGeomMap<SpatialDomains::PrismGeom>())
            {
                // Check face rotation
                if (pri->GetFace(1)->GetVid(2) != pri->GetVid(4))
                {
                    std::cerr
                        << "ERROR: Face " << pri->GetFid(1) << " (vert "
                        << pri->GetFace(1)->GetVid(2)
                        << ") not aligned to face 1 singular vert of Prism "
                        << pri->GetGlobalID() << " (vert " << pri->GetVid(4)
                        << ")" << std::endl;
                    NoRotateIssues = false;
                }
 
                // Check face rotation
                if (pri->GetFace(3)->GetVid(2) != pri->GetVid(5))
                {
                    std::cerr
                        << "ERROR: Face " << pri->GetFid(3) << " (vert "
                        << pri->GetFace(3)->GetVid(2)
                        << ") not aligned to face 3 singular vert of Prism "
                        << pri->GetGlobalID() << " (vert " << pri->GetVid(5)
                        << ")" << std::endl;
                    NoRotateIssues = false;
                }
            }
 
            if (NoRotateIssues)
            {
                std::cout << "All Prism Tri faces are correctly aligned"
                          << std::endl;
            }
 
            for ([[maybe_unused]] auto [id, hex] :
                 mesh->GetGeomMap<SpatialDomains::HexGeom>())
            {
                // Put hex checks in here
            }
        }
 
        break;
        default:
            ASSERTL0(false, "Expansion dimension not recognised");
            break;
    }
 
    //-----------------------------------------------
 
    return 0;
}

References ASSERTL0, CheckTetRotation(), Nektar::LibUtilities::SessionReader::CreateInstance(), and Nektar::SpatialDomains::MeshGraphIO::Read().