diff --git a/.gitignore b/.gitignore index 378eac25d..916fce7b8 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ build +settings.json \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 207ca93f2..21880f4e1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,10 @@ if(USE_XSDK_DEFAULTS) xsdk_compiler_flags() endif() +# require c++14 +option(ENABLE_CGNS "Enable the CGNS reader: requires c++14 extensions" OFF) +message(STATUS "ENABLE_CGNS: ${ENABLE_CGNS}") + # Set some default compiler flags that should always be used if(NOT USE_XSDK_DEFAULTS) bob_set_shared_libs() @@ -29,7 +33,15 @@ if(NOT USE_XSDK_DEFAULTS) bob_end_cxx_flags() set(CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS}") if(SCOREC_ENABLE_CXX11) - bob_cxx11_flags() + if(ENABLE_CGNS) + bob_cxx14_flags() + else() + bob_cxx14_flags() + endif() + else() + if(ENABLE_CGNS) + bob_cxx14_flags() + endif() endif() endif() message(STATUS "CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}") @@ -91,6 +103,7 @@ option(ENABLE_SIMMETRIX "Build with Simmetrix support" OFF) message(STATUS "ENABLE_SIMMETRIX: ${ENABLE_SIMMETRIX}") option(ENABLE_OMEGA_H "Enable the Omega_h interface" OFF) message(STATUS "ENABLE_OMEGA_H: ${ENABLE_OMEGA_H}") + if(ENABLE_SIMMETRIX) add_definitions(-DHAVE_SIMMETRIX) endif() @@ -117,6 +130,15 @@ if(ENABLE_OMEGA_H) bob_public_dep(Omega_h) endif() +if(ENABLE_CGNS) + find_package(CGNS) + include_directories(SYSTEM ${CGNS_INCLUDE_DIR}) + # + find_package(HDF5) + # + add_definitions(-DHAVE_CGNS) +endif() + # Include the SCOREC project packages add_subdirectory(lion) add_subdirectory(pcu) @@ -142,6 +164,9 @@ add_subdirectory(omega_h) # this INTERFACE target bundles all the enabled libraries together add_library(core INTERFACE) target_link_libraries(core INTERFACE ${SCOREC_EXPORTED_TARGETS}) +if(ENABLE_CGNS) + target_link_libraries(core INTERFACE ${CGNS_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) +endif() scorec_export_library(core) if(BUILD_EXES) diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index 881a9fe65..ed2bb2140 100644 --- a/apf/CMakeLists.txt +++ b/apf/CMakeLists.txt @@ -49,6 +49,10 @@ set(SOURCES apfMIS.cc ) +if(ENABLE_CGNS) + set(SOURCES ${SOURCES} apfCGNS.cc) +endif(ENABLE_CGNS) + # Package headers set(HEADERS apf.h @@ -94,6 +98,12 @@ target_link_libraries(apf mth ) + +if(ENABLE_CGNS) + message(STATUS ${CGNS_LIBRARIES}) + target_link_libraries(apf PRIVATE ${CGNS_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) +endif(ENABLE_CGNS) + scorec_export_library(apf) bob_end_subdir() diff --git a/apf/apf.cc b/apf/apf.cc index ceefdf53f..bd7259483 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -69,7 +69,7 @@ Field* makeField( FieldShape* shape, FieldData* data) { - PCU_ALWAYS_ASSERT( ! m->findField(name)); + PCU_ALWAYS_ASSERT_VERBOSE( ! m->findField(name), name); Field* f = 0; if (valueType == SCALAR) f = new ScalarField(); diff --git a/apf/apf.h b/apf/apf.h index 1ab3ee30f..d036d1a1f 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -13,6 +13,7 @@ #include "apfDynamicArray.h" #include +#include #include /** \file apf.h @@ -34,6 +35,7 @@ class Field; class Element; class Mesh; class MeshEntity; +class MeshTag; class VectorElement; /** \brief Mesh Elements represent the mesh coordinate vector field. */ typedef VectorElement MeshElement; @@ -629,6 +631,28 @@ for (t::iterator i = (w).begin(); \ for (t::const_iterator i = (w).begin(); \ (i) != (w).end(); ++(i)) +struct CGNSInfo +{ + // cgns_bc_name + std::string cgnsBCSName; + /* tag value + + Tag value holds [0, 1] as a + marker to indicate mesh_entities + within bc group. 1="in group", 0="not in group" + Tags set on vertices, edges, faces, and cells + */ + apf::MeshTag* bcsMarkerTag = nullptr; + // model dimension + int mdlDim = -1; + // model id + int mdlId = -1; +}; + +//using CGNSBCMap = std::map>>; +using CGNSBCMap = std::map>; +void writeCGNS(const char *prefix, Mesh *m, const CGNSBCMap &cgnsBCMap); + /** \brief Write a set of parallel VTK Unstructured Mesh files from an apf::Mesh * with binary (base64) encoding and zlib compression (if LION_COMPRESS=ON) * \details Nodal fields whose shape differs from the mesh shape will diff --git a/apf/apfCGNS.cc b/apf/apfCGNS.cc new file mode 100644 index 000000000..4debbb58b --- /dev/null +++ b/apf/apfCGNS.cc @@ -0,0 +1,1156 @@ +/* + * Author: Dr Andrew Parker (2019) - FGE Ltd + * + * This work is open source software, licensed under the terms of the + * BSD license as described in the LICENSE file in the top-level directory. + */ + +#include +#include "apfMesh.h" +#include "apfNumbering.h" +#include "apfNumberingClass.h" +#include "apfShape.h" +#include "apfFieldData.h" +#include +#include +// +#include +#include +#include +#include +#include +#include +#include +#include + +// === includes for safe_mkdir === +#include +#include /*required for mode_t for mkdir on some systems*/ +#include /*using POSIX mkdir call for SMB "foo/" path*/ +#include /* for checking the error from mkdir */ +// =============================== + +#ifdef HAVE_CGNS +// +#include +#include +// +#endif + +// Note: currently, even in 2D or 1D a full 3D vector and matrix are written out to the file + +namespace +{ +//https://proteustoolkit.org/capi/html/_parallel_mesh_converter_8cpp_source.html +using Count = std::pair; +static Count count(apf::Mesh *m, int dim) +{ + const int local = apf::countOwned(m, dim); + int total = local; + PCU_Add_Ints(&total, 1); // size of total array + return std::make_pair(total, local); +} + +struct CGNS +{ + int index = -1; + int base = -1; + int zone = -1; + const int phys_dim = 3; + std::string fname; +}; + +/* +* Disable this until I find out why all tags get an _tet_ or _pyr_ fieldname, makes the output file super messy, and fields seem duplicated in tags +void WriteTags(const CGNS &cgns, const std::vector> &orderedEnts, const std::vector> &ranges, const std::vector &orderedVertices, const int &vStart, const int &vEnd, apf::Mesh *m) +{ + const auto loopVertexTags = [&m](const auto &orderedEnts, const int &solIndex, const auto &inner, const auto &post, const int &start, const int &end) { + apf::DynamicArray tags; + m->getTags(tags); + for (std::size_t i = 0; i < tags.getSize(); ++i) + { + apf::MeshTag *t = tags[i]; + const int tagType = m->getTagType(t); + const int tagSize = m->getTagSize(t); + std::string tagName(m->getTagName(t)); + + if (tagName.size() > 32) + tagName.resize(32); // daft api + + // boring... replace with variant + std::vector idata; + std::vector ddata; + std::vector ldata; + + for (int ti = 0; ti < tagSize; ti++) + { + cgsize_t rmin[3]; + cgsize_t rmax[3]; + + rmin[0] = start; + rmax[0] = end; + //std::cout << start << " " << end << std::endl; + int fieldIndex = -1; + auto tagNameNew = tagName; + if (tagSize > 1) + tagNameNew += "_[" + std::to_string(ti) + "]"; + + for (const auto &e : orderedEnts) + { + if (m->hasTag(e, t) && m->isOwned(e)) + { + inner(e, t, tagSize, ti, idata, ddata, ldata, tagType); + } + } + + // ensure collectives are called by all, even if local has no data + int isize = idata.size(); + PCU_Add_Ints(&isize, 1); // size of total array + + int dsize = ddata.size(); + PCU_Add_Ints(&dsize, 1); // size of total array + + int lsize = ldata.size(); + PCU_Add_Ints(&lsize, 1); // size of total array + + // oddness of the api + rmin[1] = rmin[0]; + rmin[2] = rmin[0]; + rmax[1] = rmax[0]; + rmax[2] = rmax[0]; + + if (isize > 0 || dsize > 0 || lsize > 0) + post(solIndex, tagNameNew, idata, ddata, ldata, rmin, rmax, isize, dsize, lsize, fieldIndex); + } + } + }; + + const auto loopCellTags = [&m](const auto &orderedEnts, const int &solIndex, const auto &inner, const auto &post, const auto &ranges) { + apf::DynamicArray tags; + m->getTags(tags); + for (std::size_t i = 0; i < tags.getSize(); ++i) + { + apf::MeshTag *t = tags[i]; + const int tagType = m->getTagType(t); + const int tagSize = m->getTagSize(t); + std::string tagName(m->getTagName(t)); + + if (tagName.size() > 32) + tagName.resize(32); // daft api + + for (int ti = 0; ti < tagSize; ti++) + { + int fieldIndex = -1; + auto tagNameNew = tagName; + if (tagSize > 1) + tagNameNew += "_[" + std::to_string(ti) + "]"; + + for (std::size_t e = 0; e < orderedEnts.size(); e++) + { + // boring... replace with variant + std::vector idata; + std::vector ddata; + std::vector ldata; + + cgsize_t rmin[3]; + cgsize_t rmax[3]; + + rmin[0] = ranges[e].first; + rmax[0] = ranges[e].second; + + for (const auto &elm : orderedEnts[e]) + { + if (m->hasTag(elm, t) && m->isOwned(elm)) + { + inner(elm, t, tagSize, ti, idata, ddata, ldata, tagType); + } + } + + // ensure collectives are called by all, even if local has no data + int isize = idata.size(); + PCU_Add_Ints(&isize, 1); // size of total array + + int dsize = ddata.size(); + PCU_Add_Ints(&dsize, 1); // size of total array + + int lsize = ldata.size(); + PCU_Add_Ints(&lsize, 1); // size of total array + + // oddness of the api + rmin[1] = rmin[0]; + rmin[2] = rmin[0]; + rmax[1] = rmax[0]; + rmax[2] = rmax[0]; + + if (isize > 0 || dsize > 0 || lsize > 0) + post(solIndex, tagNameNew, idata, ddata, ldata, rmin, rmax, isize, dsize, lsize, fieldIndex); + } + } + } + }; + + const auto postLambda = [&cgns](const int &solIndex, const std::string &name, std::vector &idata, std::vector &ddata, std::vector &ldata, const cgsize_t *rmin, const cgsize_t *rmax, const int isize, const int dsize, const int lsize, int &fieldIndex) { + if (dsize > 0) + { + if (fieldIndex == -1) + { + if (cgp_field_write(cgns.index, cgns.base, cgns.zone, solIndex, CGNS_ENUMV(RealDouble), name.c_str(), &fieldIndex)) + cgp_error_exit(); + } + if (cgp_field_write_data(cgns.index, cgns.base, cgns.zone, solIndex, fieldIndex, &rmin[0], &rmax[0], + ddata.data())) + cgp_error_exit(); + } + else if (isize > 0) + { + if (fieldIndex == -1) + { + if (cgp_field_write(cgns.index, cgns.base, cgns.zone, solIndex, CGNS_ENUMV(Integer), name.c_str(), &fieldIndex)) + cgp_error_exit(); + } + + if (cgp_field_write_data(cgns.index, cgns.base, cgns.zone, solIndex, fieldIndex, &rmin[0], &rmax[0], + idata.data())) + cgp_error_exit(); + } + else if (lsize > 0) + { + if (fieldIndex == -1) + { + if (cgp_field_write(cgns.index, cgns.base, cgns.zone, solIndex, CGNS_ENUMV(LongInteger), name.c_str(), &fieldIndex)) + cgp_error_exit(); + } + if (cgp_field_write_data(cgns.index, cgns.base, cgns.zone, solIndex, fieldIndex, &rmin[0], &rmax[0], + ldata.data())) + cgp_error_exit(); + } + }; + + const auto innerLambda = [&m](apf::MeshEntity *elem, apf::MeshTag *tag, const int &tagSize, const int &ti, std::vector &idata, std::vector &ddata, std::vector &ldata, const int &tagType) { + if (tagType == apf::Mesh::TagType::DOUBLE) + { + std::vector vals(tagSize, -12345); + m->getDoubleTag(elem, tag, vals.data()); + ddata.push_back(vals[ti]); + } + else if (tagType == apf::Mesh::TagType::INT) + { + std::vector vals(tagSize, -12345); + m->getIntTag(elem, tag, vals.data()); + idata.push_back(vals[ti]); + } + else if (tagType == apf::Mesh::TagType::LONG) + { + std::vector vals(tagSize, -12345); + m->getLongTag(elem, tag, vals.data()); + ldata.push_back(vals[ti]); + } + else + { + std::cout << "Strange" << std::endl; + exit(-1); + } + }; + + int solIndex = -1; + + { + if (cg_sol_write(cgns.index, cgns.base, cgns.zone, "Vertex Tag Data", CGNS_ENUMV(Vertex), &solIndex)) + cg_error_exit(); + + loopVertexTags(orderedVertices, solIndex, innerLambda, postLambda, vStart, vEnd); + } + + { + if (cg_sol_write(cgns.index, cgns.base, cgns.zone, "Cell Tag Data", CGNS_ENUMV(CellCenter), &solIndex)) + cg_error_exit(); + + loopCellTags(orderedEnts, solIndex, innerLambda, postLambda, ranges); + } +} +*/ +void WriteFields(const CGNS &cgns, const std::vector> &orderedEnts, const std::vector> &ranges, const std::vector &orderedVertices, const int &vStart, const int &vEnd, apf::Mesh *m) +{ + const auto writeField = [&m, &cgns](apf::Field *f, const auto &orderedEnts, const int &solIndex, const auto &inner, const auto &post, const int &numComponents, const int &component, const std::string &fieldName, const int &start, const int &end, int &fieldIndex) { + std::vector data; + + cgsize_t rmin[3]; + cgsize_t rmax[3]; + + rmin[0] = start; + rmax[0] = end; + + apf::FieldDataOf *fieldData = f->getData(); + for (const auto &e : orderedEnts) + { + if (fieldData->hasEntity(e) && m->isOwned(e)) + { + inner(e, fieldData, data, numComponents, component); + } + } + + int size = data.size(); + PCU_Add_Ints(&size, 1); // size of total array + + // oddness of the api + rmin[1] = rmin[0]; + rmin[2] = rmin[0]; + rmax[1] = rmax[0]; + rmax[2] = rmax[0]; + if (size > 0) + { + if (fieldIndex == -1) + { + //std::cout << "FieldName " << fieldName << std::endl; + if (cgp_field_write(cgns.index, cgns.base, cgns.zone, solIndex, CGNS_ENUMV(RealDouble), fieldName.c_str(), &fieldIndex)) + cgp_error_exit(); + } + + post(solIndex, data, rmin, rmax, size, fieldIndex); + } + }; + + const auto loopCellFields = [&m, &writeField](const auto &orderedEnts, const int &solIndex, const auto &inner, const auto &post, const auto &ranges) { + for (int i = 0; i < m->countFields(); ++i) + { + apf::Field *f = m->getField(i); + const int numComponents = f->countComponents(); + std::string fieldName(f->getName()); + if (fieldName.size() > 32) + fieldName.resize(32); // daft api + + for (int component = 0; component < numComponents; component++) + { + int fieldIndex = -1; + auto fieldNameNew = fieldName; + if (numComponents > 1) + fieldNameNew += "_[" + std::to_string(component) + "]"; + + for (std::size_t e = 0; e < orderedEnts.size(); e++) + { + //std::cout << "CELL " << fieldNameNew << " " << fieldName << " " << component << " " << numComponents << " " << e << " " << ranges[e].first << " " << ranges[e].second << std::endl; + writeField(f, orderedEnts[e], solIndex, inner, post, numComponents, component, fieldNameNew, ranges[e].first, ranges[e].second, fieldIndex); + } + } + } + }; + + const auto loopVertexFields = [&m, &writeField](const auto &orderedEnts, const int &solIndex, const auto &inner, const auto &post, const int &vStart, const int &vEnd) { + for (int i = 0; i < m->countFields(); ++i) + { + apf::Field *f = m->getField(i); + const int numComponents = f->countComponents(); + std::string fieldName(f->getName()); + if (fieldName.size() > 32) + fieldName.resize(32); // daft api + + for (int component = 0; component < numComponents; component++) + { + int fieldIndex = -1; + auto fieldNameNew = fieldName; + if (numComponents > 1) + fieldNameNew += "_[" + std::to_string(component) + "]"; + + //std::cout << "VERTEX " << fieldNameNew << " " << fieldName << " " << component << " " << numComponents << std::endl; + writeField(f, orderedEnts, solIndex, inner, post, numComponents, component, fieldNameNew, vStart, vEnd, fieldIndex); + } + } + }; + + const auto postLambda = [&cgns](const int &solIndex, std::vector &ddata, const cgsize_t *rmin, const cgsize_t *rmax, const int &globalSize, const int &fieldIndex) { + if (globalSize > 0) + { + if (cgp_field_write_data(cgns.index, cgns.base, cgns.zone, solIndex, fieldIndex, &rmin[0], &rmax[0], + ddata.data())) + cgp_error_exit(); + } + }; + + const auto innerLambda = [](apf::MeshEntity *elem, apf::FieldDataOf *fieldData, std::vector &ddata, const int &numComponents, const int &component) { + std::vector vals(numComponents, -12345); + fieldData->get(elem, vals.data()); + //std::cout << numComponents << " " << component << " " << vals[0] << std::endl; + ddata.push_back(vals[component]); + }; + + int solIndex = -1; + + { + if (cg_sol_write(cgns.index, cgns.base, cgns.zone, "Vertex Field Data", CGNS_ENUMV(Vertex), &solIndex)) + cg_error_exit(); + + loopVertexFields(orderedVertices, solIndex, innerLambda, postLambda, vStart, vEnd); + } + + { + if (cg_sol_write(cgns.index, cgns.base, cgns.zone, "Cell Field Data", CGNS_ENUMV(CellCenter), &solIndex)) + cg_error_exit(); + + loopCellFields(orderedEnts, solIndex, innerLambda, postLambda, ranges); + } +} + +auto WriteVertices(const CGNS &cgns, apf::Mesh *m, apf::GlobalNumbering *gvn) +{ + int Cx = -1; + int Cy = -1; + int Cz = -1; + + if (cgns.phys_dim > 0) + { + if (cgp_coord_write(cgns.index, cgns.base, cgns.zone, CGNS_ENUMV(RealDouble), "CoordinateX", &Cx)) + cgp_error_exit(); + } + if (cgns.phys_dim > 1) + { + if (cgp_coord_write(cgns.index, cgns.base, cgns.zone, CGNS_ENUMV(RealDouble), "CoordinateY", &Cy)) + cgp_error_exit(); + } + if (cgns.phys_dim > 2) + { + if (cgp_coord_write(cgns.index, cgns.base, cgns.zone, CGNS_ENUMV(RealDouble), "CoordinateZ", &Cz)) + cgp_error_exit(); + } + + std::vector orderedVertices; + cgsize_t vertexMin[3]; + cgsize_t vertexMax[3]; + cgsize_t contigRange = -1; + { + std::array, 3> coords; + + vertexMin[0] = std::numeric_limits::max(); + vertexMax[0] = 0; + + apf::Vector3 point; + apf::MeshIterator *vertIter = m->begin(0); + apf::MeshEntity *vert = nullptr; + while ((vert = m->iterate(vertIter))) + { + if (m->isOwned(vert)) + { + const cgsize_t n = static_cast(apf::getNumber(gvn, vert, 0) + 1); // one based + if (contigRange == -1) + { + contigRange = n; + } + else + { + const auto predict = contigRange + 1; + + if (n != predict) + { + //std::cout << std::string("Range must be contigious for vertices ") + std::to_string(n) + " " + std::to_string(contigRange) << std::endl; + PCU_ALWAYS_ASSERT_VERBOSE(true == false, (std::string("Range must be contigious for vertices ") + std::to_string(n) + " " + std::to_string(contigRange)).c_str()); + } + else + contigRange = predict; // == n + } + + vertexMin[0] = std::min(vertexMin[0], n); + vertexMax[0] = std::max(vertexMax[0], n); + + m->getPoint(vert, 0, point); + coords[0].push_back(point[0]); + coords[1].push_back(point[1]); + coords[2].push_back(point[2]); + + orderedVertices.push_back(vert); + } + } + m->end(vertIter); + + // oddness of the api + vertexMin[1] = vertexMin[0]; + vertexMin[2] = vertexMin[0]; + vertexMax[1] = vertexMax[0]; + vertexMax[2] = vertexMax[0]; + + if (cgns.phys_dim > 0) + { + if (cgp_coord_write_data(cgns.index, cgns.base, cgns.zone, Cx, &vertexMin[0], &vertexMax[0], coords[0].data())) + cgp_error_exit(); + } + if (cgns.phys_dim > 1) + { + if (cgp_coord_write_data(cgns.index, cgns.base, cgns.zone, Cy, &vertexMin[0], &vertexMax[0], coords[1].data())) + cgp_error_exit(); + } + if (cgns.phys_dim > 2) + { + if (cgp_coord_write_data(cgns.index, cgns.base, cgns.zone, Cz, &vertexMin[0], &vertexMax[0], coords[2].data())) + cgp_error_exit(); + } + } + return std::make_tuple(orderedVertices, vertexMin[0], vertexMax[0]); +} + +using CellElementReturn = std::tuple>, std::vector>>; +CellElementReturn WriteElements(const CGNS &cgns, apf::Mesh *m, apf::GlobalNumbering *gvn, const int cell_dim, const Count &cellCount, const std::vector &apfElementOrder, const std::vector &cgnsElementOrder) +{ + std::vector globalNumbersByElementType(apfElementOrder.size(), 0); + std::vector numbersByElementType(apfElementOrder.size(), 0); + for (std::size_t o = 0; o < apfElementOrder.size(); o++) + { + apf::MeshIterator *cellIter = m->begin(cell_dim); + apf::MeshEntity *cell = nullptr; + int counter = 0; + while ((cell = m->iterate(cellIter))) + { + if (m->getType(cell) == apfElementOrder[o] && m->isOwned(cell)) // must be same test as below + { + counter++; + } + } + m->end(cellIter); + numbersByElementType[o] = counter; + int total = counter; + PCU_Add_Ints(&total, 1); // size of total array + globalNumbersByElementType[o] = total; + } + cgsize_t allTotal = std::accumulate(globalNumbersByElementType.begin(), globalNumbersByElementType.end(), 0); + PCU_ALWAYS_ASSERT_VERBOSE(allTotal == cellCount.first, ("Must be equal " + std::to_string(allTotal) + " " + std::to_string(cellCount.first)).c_str()); + + int globalStart = 1; // one-based + std::vector> orderedElements(apfElementOrder.size()); + std::vector> ranges(apfElementOrder.size()); + for (std::size_t o = 0; o < apfElementOrder.size(); o++) + { + std::vector elementVertices; + apf::MeshIterator *cellIter = m->begin(cell_dim); + apf::MeshEntity *cell = nullptr; + apf::Downward verts; + + while ((cell = m->iterate(cellIter))) + { + if (m->getType(cell) == apfElementOrder[o] && m->isOwned(cell)) // must be same test as above + { + const auto numVerts = m->getDownward(cell, 0, verts); + for (int i = 0; i < numVerts; i++) + { + const auto n = apf::getNumber(gvn, verts[i], 0); + elementVertices.push_back(n + 1); // one-based + } + orderedElements[o].push_back(cell); + } + } + m->end(cellIter); + + if (globalNumbersByElementType[o] > 0) + { + const int globalEnd = globalStart + globalNumbersByElementType[o] - 1; // one-based stuff + // + int sectionNumber = -1; + const std::string name = std::string(cg_ElementTypeName(cgnsElementOrder[o])) + " " + std::to_string(globalStart) + "->" + std::to_string(globalEnd); + //std::cout << "Section Name " << name << std::endl; + if (cgp_section_write(cgns.index, cgns.base, cgns.zone, name.c_str(), cgnsElementOrder[o], globalStart, globalEnd, 0, §ionNumber)) // global start, end within the file for that element type + cgp_error_exit(); + + std::vector allNumbersForThisType(PCU_Comm_Peers(), 0); + MPI_Allgather(&numbersByElementType[o], 1, MPI_INT, allNumbersForThisType.data(), 1, + MPI_INT, PCU_Get_Comm()); + + cgsize_t num = 0; + for (int i = 0; i < PCU_Comm_Self(); i++) + num += allNumbersForThisType[i]; + + cgsize_t elStart = globalStart + num; + cgsize_t elEnd = elStart + numbersByElementType[o] - 1; // one-based stuff + if (cgp_elements_write_data(cgns.index, cgns.base, cgns.zone, sectionNumber, elStart, elEnd, // per processor within the range[start, end] + elementVertices.data())) + cgp_error_exit(); + + //std::cout << std::flush << "RANK: " << PCU_Comm_Self() << " ==> " << globalStart << " " << globalEnd << " elStart " << elStart << " elEnd " << elEnd << " numbersByElementType[o] " << numbersByElementType[o] << " of " << std::string(cg_ElementTypeName(cgnsElementOrder[o])) << std::flush << std::endl; + + globalStart += globalNumbersByElementType[o]; + + ranges[o] = std::make_pair(elStart, elEnd); + } + else + ranges[o] = std::make_pair(-1, -1); + } + return std::make_tuple(orderedElements, ranges); +} + +void AddBocosToMainBase(const CGNS &cgns, const CellElementReturn &cellResults, const int &cellCount, apf::Mesh *m, const apf::CGNSBCMap &cgnsBCMap, const std::map &apf2cgns, apf::GlobalNumbering *gvn) +{ + const auto EdgeLoop = [&m](const auto &lambda, apf::MeshTag *edgeTag) { + apf::MeshIterator *edgeIter = m->begin(1); + apf::MeshEntity *edge = nullptr; + int vals[1]; + vals[0] = -1; + + while ((edge = m->iterate(edgeIter))) + { + m->getIntTag(edge, edgeTag, vals); + if (vals[0] == 1 && m->isOwned(edge)) + { + lambda(edge); + } + } + m->end(edgeIter); + }; + + const auto FaceLoop = [&m](const auto &lambda, apf::MeshTag *faceTag) { + apf::MeshIterator *faceIter = m->begin(2); + apf::MeshEntity *face = nullptr; + int vals[1]; + vals[0] = -1; + + while ((face = m->iterate(faceIter))) + { + m->getIntTag(face, faceTag, vals); + if (vals[0] == 1 && m->isOwned(face)) + { + lambda(face); + } + } + m->end(faceIter); + }; + + const auto BCEntityAdder = [&apf2cgns, &m, &cgns, &gvn](const auto &Looper, const auto &bcGroup, int &startingLocation) { + std::map> bcEntTypes; + for (const auto &r : apf2cgns) + bcEntTypes.insert(std::make_pair(r.first, std::vector())); + + const auto lambda = [&m, &bcEntTypes](apf::MeshEntity *entity) { + const auto type = m->getType(entity); + auto iter = bcEntTypes.find(type); + if (iter != bcEntTypes.end()) + { + iter->second.push_back(entity); + } + else + { + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "must not come in here"); + } + }; + // + Looper(lambda, bcGroup.bcsMarkerTag); + // + const int cacheStart = startingLocation + 1; + int cacheEnd = -1; + for (const auto &bc : bcEntTypes) + { + int startOfBCBlock = startingLocation + 1; + const int number = bc.second.size(); + int total = number; + PCU_Add_Ints(&total, 1); // size of total array + if (total > 0) + { + const auto allEnd = startOfBCBlock + total - 1; //one-based + int sectionNumber = -1; + const std::string name = bcGroup.cgnsBCSName + " " + std::to_string(startOfBCBlock) + "->" + std::to_string(allEnd); //cg_ElementTypeName(apf2cgns[bc.first]); + if (cgp_section_write(cgns.index, cgns.base, cgns.zone, name.c_str(), apf2cgns.at(std::get<0>(bc)), startOfBCBlock, allEnd, 0, + §ionNumber)) + cgp_error_exit(); + + std::vector elementVertices; + for (std::size_t e = 0; e < bc.second.size(); e++) + { + apf::Downward verts; + const auto numVerts = m->getDownward(bc.second[e], 0, verts); + for (int i = 0; i < numVerts; i++) + { + const auto n = apf::getNumber(gvn, verts[i], 0); + elementVertices.push_back(n + 1); // one-based + } + } + + std::vector allNumbersForThisType(PCU_Comm_Peers(), 0); + MPI_Allgather(&number, 1, MPI_INT, allNumbersForThisType.data(), 1, + MPI_INT, PCU_Get_Comm()); + + cgsize_t num = 0; + for (int i = 0; i < PCU_Comm_Self(); i++) + num += allNumbersForThisType[i]; + + cgsize_t elStart = startOfBCBlock + num; + cgsize_t elEnd = elStart + number - 1; // one-based stuff + + if (number == 0) + { + if (cgp_elements_write_data(cgns.index, cgns.base, cgns.zone, sectionNumber, elStart, elEnd, nullptr)) + cgp_error_exit(); + } + else + { + if (cgp_elements_write_data(cgns.index, cgns.base, cgns.zone, sectionNumber, elStart, elEnd, + elementVertices.data())) + cgp_error_exit(); + } + + // Not parallel correct + startingLocation = allEnd; + cacheEnd = allEnd; + } + } + std::vector cacheStarts(PCU_Comm_Peers(), 0); + MPI_Allgather(&cacheStart, 1, MPI_INT, cacheStarts.data(), 1, + MPI_INT, PCU_Get_Comm()); + std::vector cacheEnds(PCU_Comm_Peers(), 0); + MPI_Allgather(&cacheEnd, 1, MPI_INT, cacheEnds.data(), 1, + MPI_INT, PCU_Get_Comm()); + return std::make_pair(cacheStarts, cacheEnds); + }; + + const auto globalElementList = [](const std::vector &bcList, std::vector &allElements) { + std::vector sizes(PCU_Comm_Peers(), 0); // important initialiser + const int l = bcList.size(); + MPI_Allgather(&l, 1, MPI_INT, sizes.data(), 1, + MPI_INT, PCU_Get_Comm()); + + int totalLength = 0; + for (const auto &i : sizes) + totalLength += i; + + std::vector displacement(PCU_Comm_Peers(), -1); + displacement[0] = 0; + for (std::size_t i = 1; i < displacement.size(); i++) + displacement[i] = displacement[i - 1] + sizes[i - 1]; + + allElements.resize(totalLength); + MPI_Allgatherv(bcList.data(), bcList.size(), MPI_INT, allElements.data(), + sizes.data(), displacement.data(), MPI_INT, + PCU_Get_Comm()); + }; + + const auto doVertexBC = [&](const auto &iter) { + for (const auto &p : iter->second) + { + std::vector bcList; + apf::MeshIterator *vertIter = m->begin(0); + apf::MeshEntity *vert = nullptr; + int vals[1]; + vals[0] = -1; + + while ((vert = m->iterate(vertIter))) + { + m->getIntTag(vert, p.bcsMarkerTag, vals); + if (vals[0] == 1 && m->isOwned(vert)) + { + const auto n = apf::getNumber(gvn, vert, 0); + bcList.push_back(n + 1); // one-based + } + } + m->end(vertIter); + + // Annoying that the API requires this! + std::vector allElements; + globalElementList(bcList, allElements); + + int bcIndex = -1; + if (cg_boco_write(cgns.index, cgns.base, cgns.zone, p.cgnsBCSName.c_str(), CGNS_ENUMV(BCGeneral), CGNS_ENUMV(PointList), allElements.size(), + allElements.data(), &bcIndex)) + cg_error_exit(); + + CGNS_ENUMT(GridLocation_t) + location = CGNS_ENUMV(Vertex); + if (cg_boco_gridlocation_write(cgns.index, cgns.base, cgns.zone, bcIndex, location)) + cg_error_exit(); + } + }; + + const auto doEdgeBC = [&](const auto &iter, int &startingLocation) { + for (const auto &p : iter->second) + { + const auto se = BCEntityAdder(EdgeLoop, p, startingLocation); + for (int i = 0; i < PCU_Comm_Peers(); i++) + { + PCU_ALWAYS_ASSERT_VERBOSE(se.first[i] == se.first[0], "Must all be the same "); + PCU_ALWAYS_ASSERT_VERBOSE(se.second[i] == se.second[0], "Must all be the same "); + } + + const std::array bcRange = {{se.first[0], se.second[0]}}; + int bcIndex = -1; + if (cg_boco_write(cgns.index, cgns.base, cgns.zone, p.cgnsBCSName.c_str(), CGNS_ENUMV(BCGeneral), CGNS_ENUMV(PointRange), 2, + bcRange.data(), &bcIndex)) + cg_error_exit(); + + CGNS_ENUMT(GridLocation_t) + location = CGNS_ENUMV(EdgeCenter); + if (cg_boco_gridlocation_write(cgns.index, cgns.base, cgns.zone, bcIndex, location)) + cg_error_exit(); + } + }; + + const auto doFaceBC = [&](const auto &iter, int &startingLocation) { + for (const auto &p : iter->second) + { + const auto se = BCEntityAdder(FaceLoop, p, startingLocation); + + for (int i = 0; i < PCU_Comm_Peers(); i++) + { + PCU_ALWAYS_ASSERT_VERBOSE(se.first[i] == se.first[0], "Must all be the same "); + PCU_ALWAYS_ASSERT_VERBOSE(se.second[i] == se.second[0], "Must all be the same "); + } + const std::array bcRange = {{se.first[0], se.second[0]}}; + int bcIndex = -1; + if (cg_boco_write(cgns.index, cgns.base, cgns.zone, p.cgnsBCSName.c_str(), CGNS_ENUMV(BCGeneral), CGNS_ENUMV(PointRange), 2, + bcRange.data(), &bcIndex)) + cg_error_exit(); + + CGNS_ENUMT(GridLocation_t) + location = CGNS_ENUMV(FaceCenter); + if (cg_boco_gridlocation_write(cgns.index, cgns.base, cgns.zone, bcIndex, location)) + cg_error_exit(); + } + }; + + const auto doCellBC = [&](const auto &iter, const int &) { + for (const auto &p : iter->second) + { + std::vector bcList; + int vals[1]; + vals[0] = -1; + + const auto &elements = std::get<0>(cellResults); + const auto &ranges = std::get<1>(cellResults); + + for (std::size_t e = 0; e < elements.size(); e++) + { + std::size_t counter = ranges[e].first; + for (const auto &elm : elements[e]) + { + m->getIntTag(elm, p.bcsMarkerTag, vals); + if (vals[0] == 1 && m->isOwned(elm)) + { + bcList.push_back(counter); + } + counter++; + } + } + + // Annoying that the API requires this! + std::vector allElements; + globalElementList(bcList, allElements); + + int bcIndex = -1; + if (cg_boco_write(cgns.index, cgns.base, cgns.zone, p.cgnsBCSName.c_str(), CGNS_ENUMV(BCGeneral), CGNS_ENUMV(PointList), allElements.size(), + allElements.data(), &bcIndex)) + cg_error_exit(); + + CGNS_ENUMT(GridLocation_t) + location = CGNS_ENUMV(CellCenter); + if (cg_boco_gridlocation_write(cgns.index, cgns.base, cgns.zone, bcIndex, location)) + cg_error_exit(); + } + }; + + const auto cell_dim = m->getDimension(); + if (cell_dim == 3) + { + int startingLocation = cellCount; + + auto iter = cgnsBCMap.find("Vertex"); + if (iter != cgnsBCMap.end()) + { + doVertexBC(iter); + } + iter = cgnsBCMap.find("EdgeCenter"); + if (iter != cgnsBCMap.end()) + { + doEdgeBC(iter, startingLocation); + } + iter = cgnsBCMap.find("FaceCenter"); + if (iter != cgnsBCMap.end()) + { + doFaceBC(iter, startingLocation); + } + iter = cgnsBCMap.find("CellCenter"); + if (iter != cgnsBCMap.end()) + { + doCellBC(iter, 3); + } + } + else if (cell_dim == 2) + { + int startingLocation = cellCount; + + auto iter = cgnsBCMap.find("Vertex"); + if (iter != cgnsBCMap.end()) + { + doVertexBC(iter); + } + iter = cgnsBCMap.find("EdgeCenter"); + if (iter != cgnsBCMap.end()) + { + doEdgeBC(iter, startingLocation); + } + iter = cgnsBCMap.find("CellCenter"); + if (iter != cgnsBCMap.end()) + { + doCellBC(iter, 2); + } + } + else if (cell_dim == 1) + { + auto iter = cgnsBCMap.find("Vertex"); + if (iter != cgnsBCMap.end()) + { + doVertexBC(iter); + } + iter = cgnsBCMap.find("CellCenter"); + if (iter != cgnsBCMap.end()) + { + doCellBC(iter, 1); + } + } +} + +// copy is deliberate +void Write3DFaces(CGNS cgns, apf::Mesh *m, const Count &faceCount, const Count &vertexCount, apf::GlobalNumbering *gvn) +{ + std::array sizes; + sizes[0] = vertexCount.first; // global + sizes[1] = faceCount.first; // global + sizes[2] = 0; // nodes are unsorted, as defined by api + + { + std::string baseName("FaceMeshBase"); + if (cg_base_write(cgns.index, baseName.c_str(), 2, cgns.phys_dim, &cgns.base)) + cg_error_exit(); + } + // Write the default units at the top. + if (cg_goto(cgns.index, cgns.base, "end")) + cg_error_exit(); + + if (cg_units_write(CGNS_ENUMV(Kilogram), CGNS_ENUMV(Meter), CGNS_ENUMV(Second), CGNS_ENUMV(Kelvin), + CGNS_ENUMV(Degree))) + cg_error_exit(); + + if (cg_dataclass_write(CGNS_ENUMV(Dimensional))) + cg_error_exit(); + + { + std::string zoneName("FaceMeshZone"); + if (cg_zone_write(cgns.index, cgns.base, zoneName.c_str(), sizes.data(), CGNS_ENUMV(Unstructured), &cgns.zone)) + cg_error_exit(); + } + + // Stupid api why in the name would you have to repeat ALL vertices again.... + auto &&vertResult = WriteVertices(cgns, m, gvn); + // + std::vector apfElementOrder; + apfElementOrder.push_back(apf::Mesh::QUAD); + apfElementOrder.push_back(apf::Mesh::TRIANGLE); + + std::vector cgnsElementOrder; + cgnsElementOrder.push_back(CGNS_ENUMV(QUAD_4)); + cgnsElementOrder.push_back(CGNS_ENUMV(TRI_3)); + // + auto &&cellResult = WriteElements(cgns, m, gvn, 2, faceCount, apfElementOrder, cgnsElementOrder); + + apf::GlobalNumbering *gcn = nullptr; + gcn = apf::makeGlobal( + apf::numberOwnedDimension(m, "2D element-nums", 2)); + synchronize(gcn); + // + //WriteTags(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + // + WriteFields(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + + destroyGlobalNumbering(gcn); +} + +// copy is deliberate +void Write2DEdges(CGNS cgns, apf::Mesh *m, const Count &edgeCount, const Count &vertexCount, apf::GlobalNumbering *gvn) +{ + std::array sizes; + sizes[0] = vertexCount.first; // global + sizes[1] = edgeCount.first; // global + sizes[2] = 0; // nodes are unsorted, as defined by api + + { + std::string baseName("EdgeMeshBase"); + if (cg_base_write(cgns.index, baseName.c_str(), 1, cgns.phys_dim, &cgns.base)) + cg_error_exit(); + } + // Write the default units at the top. + if (cg_goto(cgns.index, cgns.base, "end")) + cg_error_exit(); + + if (cg_units_write(CGNS_ENUMV(Kilogram), CGNS_ENUMV(Meter), CGNS_ENUMV(Second), CGNS_ENUMV(Kelvin), + CGNS_ENUMV(Degree))) + cg_error_exit(); + + if (cg_dataclass_write(CGNS_ENUMV(Dimensional))) + cg_error_exit(); + + { + std::string zoneName("EdgeMeshZone"); + if (cg_zone_write(cgns.index, cgns.base, zoneName.c_str(), sizes.data(), CGNS_ENUMV(Unstructured), &cgns.zone)) + cg_error_exit(); + } + + // Stupid api why in the name would you have to repeat ALL vertices again.... + auto &&vertResult = WriteVertices(cgns, m, gvn); + // + std::vector apfElementOrder; + apfElementOrder.push_back(apf::Mesh::EDGE); + + std::vector cgnsElementOrder; + cgnsElementOrder.push_back(CGNS_ENUMV(BAR_2)); + // + auto &&cellResult = WriteElements(cgns, m, gvn, 1, edgeCount, apfElementOrder, cgnsElementOrder); + + apf::GlobalNumbering *gcn = nullptr; + gcn = apf::makeGlobal( + apf::numberOwnedDimension(m, "1D element-nums", 1)); + synchronize(gcn); + // + //WriteTags(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + // + WriteFields(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + + destroyGlobalNumbering(gcn); +} + +// Todo split this out into a list of calls to local functions to show process/work flow +void WriteCGNS(const char *prefix, apf::Mesh *m, const apf::CGNSBCMap &cgnsBCMap) +{ + static_assert(std::is_same::value, "cgsize_t not compiled as int"); + + const auto myRank = PCU_Comm_Self(); + const Count vertexCount = count(m, 0); + const Count edgeCount = count(m, 1); + const Count faceCount = count(m, 2); + const auto cell_dim = m->getDimension(); + const Count cellCount = count(m, cell_dim); + // for (int i = 0; i < PCU_Comm_Peers(); ++i) + // { + // if (i == PCU_Comm_Self()) + // { + // std::cout << "*******Local Mesh Stats******************\n"; + // std::cout << "Rank: " << myRank << ": Number of cells " << cellCount.second << "\n"; + // std::cout << "Rank: " << myRank << ": Number of vertices " << vertexCount.second << "\n"; + // std::cout << "Rank: " << myRank << ": Number of faces " << faceCount.second << "\n"; + // std::cout << "Rank: " << myRank << ": Number of edges " << edgeCount.second << "\n"; + // std::cout << "*****************************************\n"; + // } + // PCU_Barrier(); + // } + + PCU_Barrier(); + if (myRank == 0) + { + std::cout << "*******Global Mesh Stats*****************\n"; + std::cout << ": Number of cells " << cellCount.first << "\n"; + std::cout << ": Number of vertices " << vertexCount.first << "\n"; + std::cout << ": Number of faces " << faceCount.first << "\n"; + std::cout << ": Number of edges " << edgeCount.first << "\n"; + std::cout << "*****************************************\n"; + } + + std::array sizes; + sizes[0] = vertexCount.first; // global + sizes[1] = cellCount.first; // global + sizes[2] = 0; // nodes are unsorted, as defined by api + + // Copy communicator + auto communicator = PCU_Get_Comm(); + cgp_mpi_comm(communicator); + // + cgp_pio_mode(CGNS_ENUMV(CGP_INDEPENDENT)); + + CGNS cgns; + cgns.fname = std::string(prefix); + if (cgp_open(prefix, CGNS_ENUMV(CG_MODE_WRITE), &cgns.index)) + cgp_error_exit(); + + { + std::string baseName("VolumeMeshBase"); + if (cg_base_write(cgns.index, baseName.c_str(), cell_dim, cgns.phys_dim, &cgns.base)) + cg_error_exit(); + } + // Write the default units at the top. + if (cg_goto(cgns.index, cgns.base, "end")) + cg_error_exit(); + + if (cg_units_write(CGNS_ENUMV(Kilogram), CGNS_ENUMV(Meter), CGNS_ENUMV(Second), CGNS_ENUMV(Kelvin), + CGNS_ENUMV(Degree))) + cg_error_exit(); + + if (cg_dataclass_write(CGNS_ENUMV(Dimensional))) + cg_error_exit(); + + { + std::string zoneName("VolumeMeshZone"); + if (cg_zone_write(cgns.index, cgns.base, zoneName.c_str(), sizes.data(), CGNS_ENUMV(Unstructured), &cgns.zone)) + cg_error_exit(); + } + + apf::GlobalNumbering *gvn = nullptr; + gvn = apf::makeGlobal(apf::numberOwnedNodes(m, "node-nums")); + apf::synchronize(gvn); + // + auto &&vertResult = WriteVertices(cgns, m, gvn); + // + std::vector apfElementOrder; + apfElementOrder.push_back(apf::Mesh::HEX); + apfElementOrder.push_back(apf::Mesh::TET); + apfElementOrder.push_back(apf::Mesh::PYRAMID); + apfElementOrder.push_back(apf::Mesh::QUAD); + apfElementOrder.push_back(apf::Mesh::TRIANGLE); + apfElementOrder.push_back(apf::Mesh::EDGE); + + std::vector cgnsElementOrder; + cgnsElementOrder.push_back(CGNS_ENUMV(HEXA_8)); + cgnsElementOrder.push_back(CGNS_ENUMV(TETRA_4)); + cgnsElementOrder.push_back(CGNS_ENUMV(PYRA_5)); + cgnsElementOrder.push_back(CGNS_ENUMV(QUAD_4)); + cgnsElementOrder.push_back(CGNS_ENUMV(TRI_3)); + cgnsElementOrder.push_back(CGNS_ENUMV(BAR_2)); + // + auto &&cellResult = WriteElements(cgns, m, gvn, m->getDimension(), cellCount, apfElementOrder, cgnsElementOrder); + // + apf::GlobalNumbering *gcn = nullptr; + gcn = apf::makeGlobal(apf::numberElements(m, "3D element-nums")); + // Could be a boco at any dim + std::map apf2cgns; + apf2cgns.insert(std::make_pair(apf::Mesh::HEX, CGNS_ENUMV(HEXA_8))); + apf2cgns.insert(std::make_pair(apf::Mesh::TET, CGNS_ENUMV(TETRA_4))); + apf2cgns.insert(std::make_pair(apf::Mesh::PYRAMID, CGNS_ENUMV(PYRA_5))); + apf2cgns.insert(std::make_pair(apf::Mesh::QUAD, CGNS_ENUMV(QUAD_4))); + apf2cgns.insert(std::make_pair(apf::Mesh::TRIANGLE, CGNS_ENUMV(TRI_3))); + apf2cgns.insert(std::make_pair(apf::Mesh::EDGE, CGNS_ENUMV(BAR_2))); + // + AddBocosToMainBase(cgns, cellResult, cellCount.first, m, cgnsBCMap, apf2cgns, gvn); + // + //WriteTags(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + // + WriteFields(cgns, std::get<0>(cellResult), std::get<1>(cellResult), std::get<0>(vertResult), std::get<1>(vertResult), std::get<2>(vertResult), m); + // + if (cell_dim == 3) + { + Write3DFaces(cgns, m, faceCount, vertexCount, gvn); + // + Write2DEdges(cgns, m, edgeCount, vertexCount, gvn); + } + else if (cell_dim == 2) + { + Write2DEdges(cgns, m, edgeCount, vertexCount, gvn); + } + // + destroyGlobalNumbering(gvn); + destroyGlobalNumbering(gcn); + // + cgp_close(cgns.index); +} +} // namespace + +namespace apf +{ + +void writeCGNS(const char *prefix, Mesh *m, const apf::CGNSBCMap &cgnsBCMap) +{ +#ifdef HAVE_CGNS + WriteCGNS(prefix, m, cgnsBCMap); +#else + PCU_ALWAYS_ASSERT_VERBOSE(true == false, + "Build with ENABLE_CGNS to allow this functionality."); + exit(EXIT_FAILURE); +#endif +} + +} // namespace apf diff --git a/apf/apfConstruct.cc b/apf/apfConstruct.cc index f0f3f9b65..8a96b2b00 100644 --- a/apf/apfConstruct.cc +++ b/apf/apfConstruct.cc @@ -20,19 +20,21 @@ static void constructVerts( result[conn[i]] = m->createVert_(interior); } -static void constructElements( +static NewElements constructElements( Mesh2* m, const Gid* conn, int nelem, int etype, GlobalToVert& globalToVert) { ModelEntity* interior = m->findModelEntity(m->getDimension(), 0); int nev = apf::Mesh::adjacentCount[etype][0]; + NewElements newElements; for (int i = 0; i < nelem; ++i) { Downward verts; int offset = i * nev; for (int j = 0; j < nev; ++j) verts[j] = globalToVert[conn[j + offset]]; - buildElement(m, interior, etype, verts); + newElements.push_back(buildElement(m, interior, etype, verts)); } + return newElements; } static Gid getMax(const GlobalToVert& globalToVert) @@ -144,11 +146,11 @@ static void constructRemotes(Mesh2* m, GlobalToVert& globalToVert) } } -void assemble(Mesh2* m, const int* conn, int nelem, int etype, +NewElements assemble(Mesh2* m, const int* conn, int nelem, int etype, GlobalToVert& globalToVert) { constructVerts(m, conn, nelem, etype, globalToVert); - constructElements(m, conn, nelem, etype, globalToVert); + return constructElements(m, conn, nelem, etype, globalToVert); } void finalise(Mesh2* m, GlobalToVert& globalToVert) @@ -159,11 +161,12 @@ void finalise(Mesh2* m, GlobalToVert& globalToVert) m->acceptChanges(); } -void construct(Mesh2* m, const int* conn, int nelem, int etype, +NewElements construct(Mesh2* m, const int* conn, int nelem, int etype, GlobalToVert& globalToVert) { - assemble(m, conn, nelem, etype, globalToVert); + const auto& newElements = assemble(m, conn, nelem, etype, globalToVert); finalise(m, globalToVert); + return newElements; } void setCoords(Mesh2* m, const double* coords, int nverts, diff --git a/apf/apfConvert.h b/apf/apfConvert.h index e5209522b..fcee947f8 100644 --- a/apf/apfConvert.h +++ b/apf/apfConvert.h @@ -12,6 +12,7 @@ \brief algorithms for mesh format conversion */ #include +#include namespace apf { @@ -19,6 +20,7 @@ class Mesh; class Mesh2; class ModelEntity; class MeshEntity; +using NewElements = std::vector; /** \brief convert one mesh data structure to another \details this function will fill in a structure that fully @@ -36,7 +38,7 @@ typedef std::map GlobalToVert; assemble and finalise. The premise of assemble being that it is called multiple times for a given cell type, across several different cell types in the input mesh. */ -void assemble(Mesh2* m, const int* conn, int nelem, int etype, +NewElements assemble(Mesh2* m, const int* conn, int nelem, int etype, GlobalToVert& globalToVert); /** \brief finalise construction of a mixed-cell-type mesh from just a connectivity array @@ -63,7 +65,7 @@ void finalise(Mesh2* m, GlobalToVert& globalToVert); Note that all vertices will have zero coordinates, so it is often good to use apf::setCoords after this. */ -void construct(Mesh2* m, const int* conn, int nelem, int etype, +NewElements construct(Mesh2* m, const int* conn, int nelem, int etype, GlobalToVert& globalToVert); /** \brief Assign coordinates to the mesh diff --git a/apf/pkg_tribits.cmake b/apf/pkg_tribits.cmake index a0c6cff0c..9e2d3d048 100644 --- a/apf/pkg_tribits.cmake +++ b/apf/pkg_tribits.cmake @@ -51,6 +51,10 @@ set(APF_SOURCES apfFile.cc ) +if(ENABLE_CGNS) + set(APF_SOURCES ${APF_SOURCES} apfCGNS.cc) +endif(ENABLE_CGNS) + set(APF_HEADERS apf.h apfMesh.h diff --git a/cmake/FindCGNS.cmake b/cmake/FindCGNS.cmake new file mode 100644 index 000000000..b3803b4d6 --- /dev/null +++ b/cmake/FindCGNS.cmake @@ -0,0 +1,43 @@ + + +# +# Find the native CGNS includes and library +# +# CGNS_INCLUDE_DIR - where to find cgns.h, etc. +# CGNS_LIBRARIES - List of fully qualified libraries to link against when using CGNS. +# CGNS_FOUND - Do not attempt to use CGNS if "no" or undefined. + +find_path(CGNS_INCLUDE_DIR cgnslib.h + /usr/local/include + /usr/include +) + +find_library(CGNS_LIBRARY cgns + /usr/local/lib + /usr/lib +) + +set(CGNS_FOUND "NO") +if(CGNS_INCLUDE_DIR) + if(CGNS_LIBRARY) + set( CGNS_LIBRARIES ${CGNS_LIBRARY} ) + set( CGNS_FOUND "YES" ) + endif() +endif() + +if(CGNS_FIND_REQUIRED AND NOT CGNS_FOUND) + message(SEND_ERROR "Unable to find the requested CGNS libraries.") +endif() + +# handle the QUIETLY and REQUIRED arguments and set CGNS_FOUND to TRUE if +# all listed variables are TRUE +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(CGNS DEFAULT_MSG CGNS_LIBRARY CGNS_INCLUDE_DIR) + + +mark_as_advanced( + CGNS_INCLUDE_DIR + CGNS_LIBRARY +) + + diff --git a/cmake/bob.cmake b/cmake/bob.cmake index 90be9078b..74c5400fd 100644 --- a/cmake/bob.cmake +++ b/cmake/bob.cmake @@ -86,6 +86,18 @@ function(bob_cxx11_flags) set(CMAKE_CXX_FLAGS "${FLAGS}" PARENT_SCOPE) endfunction(bob_cxx11_flags) +function(bob_cxx14_flags) + set(FLAGS "${CMAKE_CXX_FLAGS}") + # clang only: -Werror=return-stack-address -Werror=mismatched-tags + set(FLAGS "${FLAGS} --std=c++14 -Wall -Wextra -Wpedantic -Werror -Wno-extra-semi -Werror=unused-parameter -Wno-error=deprecated-declarations") + if(CMAKE_CXX_COMPILER_ID MATCHES "Clang") + if (${PROJECT_NAME}_CXX_WARNINGS) + set(FLAGS "${FLAGS} -Wno-c++98-compat-pedantic -Wno-c++98-compat") + endif() + endif() + set(CMAKE_CXX_FLAGS "${FLAGS}" PARENT_SCOPE) +endfunction(bob_cxx14_flags) + function(bob_end_cxx_flags) set(${PROJECT_NAME}_CXX_FLAGS "" CACHE STRING "Override all C++ compiler flags") set(${PROJECT_NAME}_EXTRA_CXX_FLAGS "" CACHE STRING "Extra C++ compiler flags") diff --git a/crv/crvQuality.cc b/crv/crvQuality.cc index 04119933a..c392a182e 100644 --- a/crv/crvQuality.cc +++ b/crv/crvQuality.cc @@ -98,7 +98,7 @@ Quality::Quality(apf::Mesh* m, int algorithm_) : PCU_ALWAYS_ASSERT(algorithm >= 0 && algorithm <= 2); order = mesh->getShape()->getOrder(); PCU_ALWAYS_ASSERT(order >= 1); -}; +} /* This work is based on the approach of Geometric Validity of high-order * lagrange finite elements, theory and practical guidance, diff --git a/mds/CMakeLists.txt b/mds/CMakeLists.txt index 785246c03..528c92525 100644 --- a/mds/CMakeLists.txt +++ b/mds/CMakeLists.txt @@ -31,6 +31,10 @@ set(SOURCES mdsUgrid.cc ) +if(ENABLE_CGNS) + set(SOURCES ${SOURCES} mdsCGNS.cc) +endif(ENABLE_CGNS) + # Package headers set(HEADERS apfMDS.h @@ -59,6 +63,11 @@ target_link_libraries(mds apf ) +if(ENABLE_CGNS) + message(STATUS ${CGNS_LIBRARIES}) + target_link_libraries(mds PRIVATE ${CGNS_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) +endif(ENABLE_CGNS) + scorec_export_library(mds) bob_end_subdir() diff --git a/mds/apfMDS.cc b/mds/apfMDS.cc index 3e8dddd94..121027d43 100644 --- a/mds/apfMDS.cc +++ b/mds/apfMDS.cc @@ -410,7 +410,7 @@ class MeshMDS : public Mesh2 MeshTag* createIntTag(const char* name, int size) { mds_tag* tag; - PCU_ALWAYS_ASSERT(!mds_find_tag(&mesh->tags, name)); + PCU_ALWAYS_ASSERT_VERBOSE(!mds_find_tag(&mesh->tags, name), name); tag = mds_create_tag(&(mesh->tags),name, sizeof(int)*size, Mesh::INT); return reinterpret_cast(tag); diff --git a/mds/apfMDS.h b/mds/apfMDS.h index ef745e576..284f85b06 100644 --- a/mds/apfMDS.h +++ b/mds/apfMDS.h @@ -31,7 +31,12 @@ \brief Interface to the compact Mesh Data Structure */ #include - +// +// AJP: including apf.h for single define: CGNSBCMap +// AJP: alternative is to allow common cgns base header +// but trying to avoid that since it's not core functionality +#include +// struct gmi_model; namespace apf { @@ -41,6 +46,7 @@ class Mesh2; class MeshTag; class MeshEntity; class Migration; +class Field; /** \brief a map from global ids to vertex objects */ typedef std::map GlobalToVert; @@ -192,6 +198,11 @@ int getMdsIndex(Mesh2* in, MeshEntity* e); so call apf::reorderMdsMesh after any mesh modification. */ MeshEntity* getMdsEntity(Mesh2* in, int dimension, int index); +Mesh2* loadMdsFromCGNS(gmi_model* g, const char* filename, CGNSBCMap& cgnsBCMap); + +// names of mesh data to read from file: (VERTEX, VelocityX; CellCentre, Pressure) +Mesh2* loadMdsFromCGNS(gmi_model* g, const char* filename, CGNSBCMap& cgnsBCMap, const std::vector>& meshData); + Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename); Mesh2* loadMdsFromUgrid(gmi_model* g, const char* filename); diff --git a/mds/mdsCGNS.cc b/mds/mdsCGNS.cc new file mode 100644 index 000000000..955e95031 --- /dev/null +++ b/mds/mdsCGNS.cc @@ -0,0 +1,1754 @@ +/* Author: Dr Andrew Parker (2019) - FGE Ltd + + Notes: these functions have been exclusively developed using meshes from Salome + +*/ + +#include "apf.h" +#include "apfMDS.h" +#include "apfMesh2.h" +#include "apfShape.h" +#include "gmi.h" /* this is for gmi_getline... */ +#include +#include +#include +#include +#include "apfFieldData.h" +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// +#ifdef HAVE_CGNS +// +#include +#include +// +#endif + +namespace +{ +#ifndef NDEBUG // debug settings, cmake double negative.... +const bool debugOutput = true; // probably will not get away with c++17 +//static constexpr bool debugOutput = true; // probably will not get away with c++17 +#else // optimised setting +const bool debugOutput = false; // probably will not get away with c++17 + //static constexpr bool debugOutput = false; // probably will not get away with c++17 +#endif + +static std::string SupportedCGNSElementTypeToString(const CGNS_ENUMT(ElementType_t) & elementType) +{ + if (elementType == CGNS_ENUMV(NODE)) + return "NODE"; + else if (elementType == CGNS_ENUMV(BAR_2)) + return "BAR_2"; + else if (elementType == CGNS_ENUMV(QUAD_4)) + return "QUAD_4"; + else if (elementType == CGNS_ENUMV(TRI_3)) + return "TRI_3"; + else if (elementType == CGNS_ENUMV(HEXA_8)) + return "HEXA_8"; + else if (elementType == CGNS_ENUMV(TETRA_4)) + return "TETRA_4"; + else if (elementType == CGNS_ENUMV(PYRA_5)) + return "PYRA_5"; + else + { + std::cout << __LINE__ << " " + << "Element type " + << " " << cg_ElementTypeName(elementType) << " not supported by apf/PUMI " << std::endl; + return ""; + } +} + +struct MeshData +{ + MeshData(const int &s, const int &f, const CGNS_ENUMT(DataType_t) & dt, const CGNS_ENUMT(GridLocation_t) & l, const std::string &n) : si(s), fi(f), datatype(dt), location(l), name(n) + { + } + + int si = -1; + int fi = -1; + CGNS_ENUMT(DataType_t) + datatype; + CGNS_ENUMT(GridLocation_t) + location; + std::string name; + bool process = true; + + bool operator==(const MeshData &other) const + { + if (si != other.si) + return false; + + if (fi != other.fi) + return false; + + if (datatype != other.datatype) + return false; + + if (location != other.location) + return false; + + if (name != other.name) + return false; + + return true; + } + + bool operator!=(const MeshData &other) const + { + return !this->operator==(other); + } +}; + +struct MeshDataGroup +{ + std::map components; + bool process = true; + + CGNS_ENUMT(GridLocation_t) + location() const + { + return components.at(0).location; + } + + CGNS_ENUMT(DataType_t) + datatype() const + { + return components.at(0).datatype; + } + + std::string name() const + { + // pattern matching component writer in apfCGNS.cc, patterns must be kept in-sync with this file + const std::string end("_["); + const std::size_t found = components.at(0).name.find(end); + if (found != std::string::npos) + return components.at(0).name.substr(0, found); + else if (components.size() == 1) + return components.at(0).name; + else + { + info(); + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "Problem with name"); + return ""; + } + } + + std::size_t size() const + { + return components.size(); + } + + int sIndex(const int index) const + { + return components.at(index).si; + } + + int fIndex(const int index) const + { + return components.at(index).fi; + } + + void insert(int i, MeshData &d) + { + components.insert(std::make_pair(i, d)); + } + + void info() const + { + if (components.size() == 1) + { + std::cout << "Scalar Group has " << components.size() << " related componenets: " << std::endl; + for (const auto m : components) + std::cout << "Field " << m.second.name << " @ " << m.second.si << " " << m.second.fi << std::endl; + } + else if (components.size() == 3) + { + std::cout << "Vector Group has " << components.size() << " related componenets: " << std::endl; + for (const auto m : components) + std::cout << "Field " << m.second.name << " @ " << m.second.si << " " << m.second.fi << std::endl; + } + else if (components.size() == 9) + { + std::cout << "Matrix Group has " << components.size() << " related componenets: " << std::endl; + for (const auto m : components) + std::cout << "Field " << m.second.name << " @ " << m.second.si << " " << m.second.fi << std::endl; + } + else + { + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "Tensor not accounted for"); + } + } +}; + +template +void DebugParallelPrinter(std::ostream &out, Arg &&arg, Args &&... args) +{ + // if constexpr (debugOutput) // probably will not get away with c++17 + if (debugOutput) + { + for (int i = 0; i < PCU_Comm_Peers(); i++) + { + if (i == PCU_Comm_Self()) + { + out << "Rank [" << i << "] " << std::forward(arg); + //((out << ", " << std::forward(args)), ...); // probably will not get away with c++17 + using expander = int[]; + (void)expander{0, (void(out << ", " << std::forward(args)), 0)...}; + out << "\n"; + out << std::flush; + } + // removed this: can't guarantee this is collectively called, order not ensured therefore. + //PCU_Barrier(); + } + } +} + +template +void Kill(const int fid, Args &&... args) +{ + DebugParallelPrinter(std::cout, "***** CGNS ERROR", args...); + + if (PCU_Comm_Initialized()) + { + cgp_close(fid); + // Finalize the MPI environment. + PCU_Comm_Free(); + MPI_Finalize(); + cgp_error_exit(); + exit(EXIT_FAILURE); + } + else + { + cg_close(fid); + cg_error_exit(); + exit(EXIT_FAILURE); + } +} + +void Kill(const int fid) +{ + DebugParallelPrinter(std::cout, "***** CGNS ERROR"); + + if (PCU_Comm_Initialized()) + { + cgp_close(fid); + // Finalize the MPI environment. + PCU_Comm_Free(); + MPI_Finalize(); + cgp_error_exit(); + exit(EXIT_FAILURE); + } + else + { + cg_close(fid); + cg_error_exit(); + exit(EXIT_FAILURE); + } +} + +auto ReadCGNSCoords(int cgid, int base, int zone, int ncoords, int nverts, const std::vector &, const apf::GlobalToVert &globalToVert) +{ + // Read min required as defined by consecutive range + // make one based as ReadElements makes zero based + const cgsize_t lowest = globalToVert.begin()->first + 1; + const cgsize_t highest = globalToVert.rbegin()->first + 1; + DebugParallelPrinter(std::cout, "From globalToVert ", lowest, highest, nverts); + + cgsize_t range_min[3]; + range_min[0] = lowest; + range_min[1] = lowest; + range_min[2] = lowest; + cgsize_t range_max[3]; + range_max[0] = highest; + range_max[1] = highest; + range_max[2] = highest; + + std::vector> ordinates(3); + const cgsize_t numToRead = range_max[0] - range_min[0] + 1; // one based + ordinates[0].resize(numToRead, 0.0); + ordinates[1].resize(numToRead, 0.0); + ordinates[2].resize(numToRead, 0.0); + + std::vector coord_names(3); + coord_names[0].resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + coord_names[1].resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + coord_names[2].resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + + for (int d = 0; d < ncoords; d++) + { + CGNS_ENUMT(DataType_t) + datatype = CGNS_ENUMV(DataTypeNull); + + if (cg_coord_info(cgid, base, zone, d + 1, &datatype, &coord_names[d][0])) + { + std::cout << __LINE__ << " CGNS is dead " << std::endl; + Kill(cgid); + } + const auto coord_name = std::string(coord_names[d].c_str()); + //boost::algorithm::trim(coord_name); // can't be bothered including boost + + auto &ordinate = ordinates[d]; + if (cgp_coord_read_data(cgid, base, zone, d + 1, &range_min[0], &range_max[0], + ordinate.data())) + { + std::cout << __LINE__ << " CGNS is dead " << std::endl; + Kill(cgid); + } + } + // to be clear, indices passed back are global, zero based + // should only return what's inside globalToVert, but above I have to read more.... + std::map> points; + int counter = lowest; + for (int i = 0; i < numToRead; i++) + { + int zeroBased = counter - 1; // remove as per the addition above + auto iter = globalToVert.find(zeroBased); + if (iter != globalToVert.end()) + { + points.insert(std::make_pair(zeroBased, std::array{{ordinates[0][i], ordinates[1][i], ordinates[2][i]}})); + } + counter++; + } + + return points; +} + +void SimpleElementPartition(std::vector &numberToReadPerProc, std::vector &startingIndex, int el_start /* one based */, int el_end, int numElements) +{ + numberToReadPerProc.resize(PCU_Comm_Peers(), 0); + const cgsize_t blockSize = static_cast( + std::floor(static_cast(numElements) / + static_cast(PCU_Comm_Peers()))); + + DebugParallelPrinter(std::cout, "BlockSize", blockSize, "numElements", numElements, "comms Size", PCU_Comm_Peers()); + + cgsize_t counter = 0; + if (blockSize == 0 && numElements > 0) + numberToReadPerProc[0] = numElements; + else + { + bool keepGoing = true; + while (keepGoing) + { + for (int p = 0; p < PCU_Comm_Peers() && keepGoing; p++) + { + if (counter + blockSize <= numElements) + { + counter += blockSize; + numberToReadPerProc[p] += blockSize; + } + else + { + const auto &delta = numElements - counter; + counter += delta; + numberToReadPerProc[p] += delta; + keepGoing = false; + } + } + } + } + DebugParallelPrinter(std::cout, "Sanity check: Counter", counter, "numElements", numElements); + + DebugParallelPrinter(std::cout, "numberToReadPerProc for rank", PCU_Comm_Self(), "is:", numberToReadPerProc[PCU_Comm_Self()]); + + startingIndex.resize(PCU_Comm_Peers(), 0); + startingIndex[0] = el_start; + for (std::size_t i = 1; i < numberToReadPerProc.size(); i++) + { + if (numberToReadPerProc[i] > 0) + startingIndex[i] = startingIndex[i - 1] + numberToReadPerProc[i - 1]; + } + + DebugParallelPrinter(std::cout, "Element start, end, numElements ", el_start, + el_end, numElements); + + DebugParallelPrinter(std::cout, "startingIndex for rank", PCU_Comm_Self(), "is:", startingIndex[PCU_Comm_Self()]); + + DebugParallelPrinter(std::cout, "Returning from SimpleElementPartition \n"); +} + +using Pair = std::pair; +using LocalElementRanges = std::vector; // one based + +auto ReadElements(int cgid, int base, int zone, int section, int el_start /* one based */, int el_end, int numElements, int verticesPerElement, LocalElementRanges &localElementRanges) +{ + std::vector numberToReadPerProc; + std::vector startingIndex; + SimpleElementPartition(numberToReadPerProc, startingIndex, el_start, el_end, numElements); + + std::vector vertexIDs; + if (numberToReadPerProc[PCU_Comm_Self()] > 0) + vertexIDs.resize(numberToReadPerProc[PCU_Comm_Self()] * verticesPerElement, + -1234567); + + const auto start = startingIndex[PCU_Comm_Self()]; + const auto end = startingIndex[PCU_Comm_Self()] + numberToReadPerProc[PCU_Comm_Self()] - 1; + DebugParallelPrinter(std::cout, "Range", start, "to", end, numberToReadPerProc[PCU_Comm_Self()]); + // + cgp_elements_read_data(cgid, base, zone, section, start, + end, vertexIDs.data()); + + if (numberToReadPerProc[PCU_Comm_Self()] > 0) + { + localElementRanges.push_back(std::make_pair(start, end)); + } + + // remove CGNS one-based offset + // to be clear these are the node ids defining the elements, not element ids + for (auto &i : vertexIDs) + { + i = i - 1; + } + + return std::make_tuple(vertexIDs, numberToReadPerProc[PCU_Comm_Self()]); +} + +struct CGNSBCMeta +{ + std::string bocoName; + CGNS_ENUMT(BCType_t) + bocoType = CGNS_ENUMV(BCTypeNull); + CGNS_ENUMT(PointSetType_t) + ptsetType = CGNS_ENUMV(PointSetTypeNull); + cgsize_t npnts = -1; + std::vector normalindices; + cgsize_t normalListSize = -1; + CGNS_ENUMT(DataType_t) + normalDataType = CGNS_ENUMV(DataTypeNull); + int ndataset = -1; + CGNS_ENUMT(GridLocation_t) + location = CGNS_ENUMV(GridLocationNull); + std::string locationName; + cgsize_t minElementId = -1; + cgsize_t maxElementId = -1; + std::vector bcElementIds; + + void Info() const + { + std::cout << "BC named: " << bocoName << ", located on: " << locationName << std::endl; + std::cout << "\tHas " << npnts << " elements on the bc stored as a " << cg_PointSetTypeName(ptsetType) << std::endl; + std::cout << "\tThe min and max Element Ids for this bcs are: " << minElementId << " " << maxElementId << std::endl; + if (debugOutput) + //if constexpr (debugOutput)// probably will not get away with c++17 + { + std::cout << "\tThe element Ids that are tagged with this bcs are: "; + for (const auto &i : bcElementIds) + std::cout << i << " "; + std::cout << std::endl; + } + } +}; + +struct BCInfo +{ + std::string bcName; // user provided + std::string cgnsLocation; // for debug + std::unordered_set vertexIds; // zero-based global to relate to GlobalToVert + apf::MeshTag *tag = nullptr; + apf::Field *field = nullptr; // debug, outputting to vtk + + void Clean(apf::Mesh2 *m) + { + m->removeField(field); + apf::destroyField(field); + apf::removeTagFromDimension(m, tag, 0); + m->destroyTag(tag); + } + + void TagVertices(const int cgid, apf::Mesh2 *m, apf::GlobalToVert &globalToVert) + { + // tm = "tag marker" for that bcName, fm = "field marker" for that bcName + tag = m->createIntTag(("debug_tm_" + bcName).c_str(), 1); // 1 is size of tag + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(0); + int vals[1]; + vals[0] = 0; + while ((elem = m->iterate(it))) + m->setIntTag(elem, tag, vals); + m->end(it); + + for (const auto &v : vertexIds) + { + auto iter = globalToVert.find(v); + if (iter != globalToVert.end()) + { + apf::MeshEntity *elem = iter->second; + vals[0] = 1; + m->setIntTag(elem, tag, vals); + } + else + { + Kill(cgid, "GlobalToVert lookup problem", v); + } + } + + if (debugOutput) + //if constexpr (debugOutput) // probably will not get away with c++17 + { + // for debug output, tags aren't written to vtk... + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(0); + field = apf::createFieldOn(m, ("debug_fm_" + bcName).c_str(), apf::SCALAR); + + int vals[1]; + double dval; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, tag, vals); + dval = vals[0]; + apf::setScalar(field, elem, 0, dval); + } + m->end(it); + } + + // Notes: + // I do not exchange the tag values (even if that can be done). + // I'm assuming in parallel all partitions that need the vertex that + // falls within a given group will mark that vertex accordingly. + // I assume therefore that vertices on a processor boundary, are marked + // by all procs that share it. + // TODO: generate test that proves this works. Done: for quads tested 4-way + } + + /** + * Meanings/concept follows: https://cgns.github.io/CGNS_docs_current/sids/misc.html section 12.4 + +---------------+--------------+--------------+--------------+-------------------------+ + | | GridLocation | GridLocation | GridLocation | GridLocation | + +---------------+--------------+--------------+--------------+-------------------------+ + | CellDimension | Vertex | EdgeCenter | *FaceCenter | CellCenter | + | 1 | vertices | - | - | cells (line elements) | + | 2 | vertices | edges | - | cells (area elements) | + | 3 | vertices | edges | faces | cells (volume elements) | + +---------------+--------------+--------------+--------------+-------------------------+ + **/ + void TagBCEntities(const int cgid, apf::Mesh2 *m, apf::CGNSBCMap &cgnsBCMap) + { + const auto Add = [&cgnsBCMap](const std::string &location, const std::string &tagName, apf::MeshTag *field) { + auto iter = cgnsBCMap.find(location); + if (iter != cgnsBCMap.end()) + { + apf::CGNSInfo info; + info.cgnsBCSName = tagName; + info.bcsMarkerTag = field; + iter->second.push_back(info); + } + else + { + std::vector infos; + apf::CGNSInfo info; + info.cgnsBCSName = tagName; + info.bcsMarkerTag = field; + infos.push_back(info); + cgnsBCMap.insert(std::make_pair(location, infos)); + } + }; + + const auto VertexLoop = [&m](const std::string &tagName, apf::MeshTag *vertexTag) { + apf::MeshTag *bcTag = nullptr; + bcTag = m->createIntTag(tagName.c_str(), 1); // 1 is size of tag + + apf::MeshIterator *vertIter = m->begin(0); + apf::MeshEntity *vert = nullptr; + int vals[1]; + vals[0] = 0; + while ((vert = m->iterate(vertIter))) + { + m->setIntTag(vert, bcTag, vals); + } + m->end(vertIter); + + vertIter = m->begin(0); + while ((vert = m->iterate(vertIter))) + { + bool allTagged = true; + m->getIntTag(vert, vertexTag, vals); + if (vals[0] == 0) + allTagged = false; + + if (allTagged) + { + vals[0] = 1; + m->setIntTag(vert, bcTag, vals); + } + } + m->end(vertIter); + + if (debugOutput) + //if constexpr (debugOutput) // probably will not get away with c++17 + { + // for debug output, tags aren't written to vtk... + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(0); + auto *field = apf::createFieldOn(m, ("debug_" + tagName).c_str(), apf::SCALAR); + + int vals[1]; + double dval; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, bcTag, vals); + dval = vals[0]; + apf::setScalar(field, elem, 0, dval); + } + m->end(it); + } + + return bcTag; + }; + + const auto EdgeLoop = [&m](const std::string &tagName, apf::MeshTag *vertexTag) { + apf::MeshTag *bcTag = nullptr; + bcTag = m->createIntTag(tagName.c_str(), 1); // 1 is size of tag + + apf::MeshIterator *edgeIter = m->begin(1); + apf::MeshEntity *edge = nullptr; + int vals[1]; + vals[0] = 0; + while ((edge = m->iterate(edgeIter))) + { + m->setIntTag(edge, bcTag, vals); + } + m->end(edgeIter); + + apf::Downward verts; + edgeIter = m->begin(1); + while ((edge = m->iterate(edgeIter))) + { + const auto numVerts = m->getDownward(edge, 0, verts); + bool allTagged = true; + for (int i = 0; i < numVerts; i++) + { + m->getIntTag(verts[i], vertexTag, vals); + if (vals[0] == 0) + allTagged = false; + } + if (allTagged) + { + vals[0] = 1; + m->setIntTag(edge, bcTag, vals); + } + } + m->end(edgeIter); + + if (debugOutput) + //if constexpr (debugOutput) // probably will not get away with c++17 + { + // for debug output, tags aren't written to vtk... + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(1); + auto *field = apf::createField(m, ("debug_" + tagName).c_str(), apf::SCALAR, apf::getConstant(1)); + + int vals[1]; + double dval; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, bcTag, vals); + dval = vals[0]; + apf::setScalar(field, elem, 0, dval); + } + m->end(it); + } + + return bcTag; + }; + + const auto FaceLoop = [&m](const std::string &tagName, apf::MeshTag *vertexTag) { + apf::MeshTag *bcTag = nullptr; + bcTag = m->createIntTag(tagName.c_str(), 1); // 1 is size of tag + + apf::MeshIterator *faceIter = m->begin(2); + apf::MeshEntity *face = nullptr; + int vals[1]; + vals[0] = 0; + while ((face = m->iterate(faceIter))) + { + m->setIntTag(face, bcTag, vals); + } + m->end(faceIter); + + apf::Downward verts; + faceIter = m->begin(2); + while ((face = m->iterate(faceIter))) + { + const auto numVerts = m->getDownward(face, 0, verts); + bool allTagged = true; + for (int i = 0; i < numVerts; i++) + { + m->getIntTag(verts[i], vertexTag, vals); + if (vals[0] == 0) + allTagged = false; + } + if (allTagged) + { + vals[0] = 1; + m->setIntTag(face, bcTag, vals); + } + } + m->end(faceIter); + + if (debugOutput) + //if constexpr (debugOutput) // probably will not get away with c++17 + { + // for debug output, tags aren't written to vtk... + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(2); + auto *field = apf::createField(m, ("debug_" + tagName).c_str(), apf::SCALAR, apf::getConstant(2)); + + int vals[1]; + double dval; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, bcTag, vals); + dval = vals[0]; + apf::setScalar(field, elem, 0, dval); + } + m->end(it); + } + + return bcTag; + }; + + const auto CellLoop = [&m](const std::string &tagName, apf::MeshTag *vertexTag, int dim) { + apf::MeshTag *bcTag = nullptr; + bcTag = m->createIntTag(tagName.c_str(), 1); // 1 is size of tag + + apf::MeshIterator *cellIter = m->begin(dim); + apf::MeshEntity *cell = nullptr; + int vals[1]; + vals[0] = 0; + while ((cell = m->iterate(cellIter))) + { + m->setIntTag(cell, bcTag, vals); + } + m->end(cellIter); + + apf::Downward verts; + cellIter = m->begin(dim); + while ((cell = m->iterate(cellIter))) + { + const auto numVerts = m->getDownward(cell, 0, verts); + bool allTagged = true; + for (int i = 0; i < numVerts; i++) + { + m->getIntTag(verts[i], vertexTag, vals); + if (vals[0] == 0) + allTagged = false; + } + if (allTagged) + { + vals[0] = 1; + m->setIntTag(cell, bcTag, vals); + } + } + m->end(cellIter); + + if (debugOutput) + //if constexpr (debugOutput) // probably will not get away with c++17 + { + // for debug output, tags aren't written to vtk... + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(dim); + auto *field = apf::createField(m, ("debug_" + tagName).c_str(), apf::SCALAR, apf::getConstant(dim)); + + int vals[1]; + double dval; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, bcTag, vals); + dval = vals[0]; + apf::setScalar(field, elem, 0, dval); + } + m->end(it); + } + + return bcTag; + }; + + if (m->getDimension() == 3) // working with a 3D mesh + { + if (cgnsLocation == "Vertex") + { + apf::MeshTag *bcTag = VertexLoop(bcName, tag); + Add("Vertex", bcName, bcTag); + } + else if (cgnsLocation == "EdgeCenter") + { + apf::MeshTag *bcTag = EdgeLoop(bcName, tag); + Add("EdgeCenter", bcName, bcTag); + } + else if (cgnsLocation == "FaceCenter") + { + apf::MeshTag *bcTag = FaceLoop(bcName, tag); + Add("FaceCenter", bcName, bcTag); + } + else if (cgnsLocation == "CellCenter") + { + apf::MeshTag *bcTag = CellLoop(bcName, tag, m->getDimension()); + Add("CellCenter", bcName, bcTag); + } + else + Kill(cgid, "Unknown BC Type", cgnsLocation); + } + else if (m->getDimension() == 2) // working with a 2D mesh + { + if (cgnsLocation == "Vertex") + { + apf::MeshTag *bcTag = VertexLoop(bcName, tag); + Add("Vertex", bcName, bcTag); + } + else if (cgnsLocation == "EdgeCenter") + { + apf::MeshTag *bcTag = EdgeLoop(bcName, tag); + Add("EdgeCenter", bcName, bcTag); + } + else if (cgnsLocation == "FaceCenter") + { + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "Can't have a FaceCenter BC in a 2D mesh"); + } + else if (cgnsLocation == "CellCenter") + { + apf::MeshTag *bcTag = CellLoop(bcName, tag, m->getDimension()); + Add("CellCenter", bcName, bcTag); + } + } + else if (m->getDimension() == 1) // working with a 1D mesh + { + if (cgnsLocation == "Vertex") + { + apf::MeshTag *bcTag = VertexLoop(bcName, tag); + Add("Vertex", bcName, bcTag); + } + else if (cgnsLocation == "EdgeCenter") + { + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "Can't have a EdgeCenter BC in a 1D mesh"); + } + else if (cgnsLocation == "FaceCenter") + { + PCU_ALWAYS_ASSERT_VERBOSE(true == false, "Can't have a FaceCenter BC in a 1D mesh"); + } + else if (cgnsLocation == "CellCenter") + { + apf::MeshTag *bcTag = CellLoop(bcName, tag, m->getDimension()); + Add("CellCenter", bcName, bcTag); + } + } + + if (!debugOutput) + Clean(m); + } +}; // namespace + +void ReadBCInfo(const int cgid, const int base, const int zone, const int nBocos, const int physDim, const int cellDim, const int nsections, std::vector &bcInfos, const apf::GlobalToVert &globalToVert) +{ + // Read the BCS. + std::vector bcMetas(nBocos); + bcInfos.resize(nBocos); + // + for (int boco = 1; boco <= nBocos; boco++) + { + auto &bcInfo = bcInfos[boco - 1]; + auto &bcMeta = bcMetas[boco - 1]; + bcMeta.normalindices.resize(physDim); + + bcMeta.bocoName.resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + bool pointRange = false; + + if (cg_boco_info(cgid, base, zone, boco, &bcMeta.bocoName[0], &bcMeta.bocoType, + &bcMeta.ptsetType, &bcMeta.npnts, bcMeta.normalindices.data(), &bcMeta.normalListSize, + &bcMeta.normalDataType, &bcMeta.ndataset)) + Kill(cgid, "Failed cg_boco_info"); + + if (bcMeta.ptsetType == CGNS_ENUMV(PointList) || (bcMeta.ptsetType == CGNS_ENUMV(PointRange))) + { + bcMeta.bocoName = std::string(bcMeta.bocoName.c_str()); + //boost::algorithm::trim(bcMeta.bocoName); // can't be bothered including boost + + if (cg_boco_gridlocation_read(cgid, base, zone, boco, &bcMeta.location)) + Kill(cgid, "Failed cg_boco_gridlocation_read"); + + bcMeta.locationName = cg_GridLocationName(bcMeta.location); + + if (bcMeta.ptsetType == CGNS_ENUMV(PointRange)) + pointRange = true; + } + else + { + Kill(cgid, + "TODO: Can only work with " + "PointList and PointRange BC " + "types at the moment"); + } + + bcMeta.minElementId = -1; + bcMeta.maxElementId = -1; + bcMeta.bcElementIds.resize(bcMeta.npnts, -1); + + // here I read say elements with bcs as: 5, 3, 7, 9, and then read below ALL elements from 3->9, + // but I only need, 3, 5, 7, 9, so don't need elements 4, 6, 8 + { + if (bcMeta.locationName == "Vertex") // && cellDim == 1) + { + if (cg_boco_read(cgid, base, zone, boco, bcMeta.bcElementIds.data(), NULL)) + Kill(cgid, "Failed cg_boco_read"); + } + else if (bcMeta.locationName == "EdgeCenter") // && cellDim == 2) + { + if (cg_boco_read(cgid, base, zone, boco, bcMeta.bcElementIds.data(), NULL)) + Kill(cgid, "Failed cg_boco_read"); + } + else if (bcMeta.locationName == "FaceCenter") // && cellDim == 3) + { + if (cg_boco_read(cgid, base, zone, boco, bcMeta.bcElementIds.data(), NULL)) + Kill(cgid, "Failed cg_boco_read"); + } + else if (bcMeta.locationName == "CellCenter") // && cellDim == 3) + { + if (cg_boco_read(cgid, base, zone, boco, bcMeta.bcElementIds.data(), NULL)) + Kill(cgid, "Failed cg_boco_read"); + } + else + Kill(cgid, "Failed Location test for BC Type", bcMeta.locationName, + cellDim); + + if (pointRange) + { + // Check this is correct, I'm just trying to fill a contiguous range from [start, end] + PCU_ALWAYS_ASSERT_VERBOSE(bcMeta.bcElementIds.size() == 2, "wrong size"); + const auto start = bcMeta.bcElementIds[0]; + const auto end = bcMeta.bcElementIds[1]; + const auto size = end - start + 1; + bcMeta.bcElementIds.resize(size, -1); + std::iota(std::begin(bcMeta.bcElementIds), std::end(bcMeta.bcElementIds), start); + } + + bcMeta.minElementId = *std::min_element(bcMeta.bcElementIds.begin(), bcMeta.bcElementIds.end()); + bcMeta.maxElementId = *std::max_element(bcMeta.bcElementIds.begin(), bcMeta.bcElementIds.end()); + bcMeta.Info(); + bcInfo.bcName = bcMeta.bocoName; + bcInfo.cgnsLocation = bcMeta.locationName; + } + + std::vector vertexIDs; + if (bcMeta.locationName != "Vertex") + { + std::unordered_set elementsToConsider; + for (const auto &j : bcMeta.bcElementIds) + elementsToConsider.insert(j); + + for (int section = 1; section <= nsections; section++) + { + std::string sectionName; + sectionName.resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + + CGNS_ENUMT(ElementType_t) + elementType = CGNS_ENUMV(ElementTypeNull); + cgsize_t el_start = -1; + cgsize_t el_end = -1; + int num_bndry = -1; + int parent_flag = -1; + //cgsize_t numElements = -1; + int verticesPerElement = -1; + + cg_section_read(cgid, base, zone, section, §ionName[0], + &elementType, &el_start, &el_end, + &num_bndry, &parent_flag); + + //numElements = el_end - el_start + 1; + + cg_npe(elementType, &verticesPerElement); + + const cgsize_t es = el_start; + const cgsize_t ee = el_end; + + cgsize_t range_min = bcMeta.minElementId; + cgsize_t range_max = bcMeta.maxElementId; + + bool doRead = true; + if (es < bcMeta.minElementId && ee < bcMeta.minElementId) + doRead = false; + if (es > bcMeta.maxElementId && es > bcMeta.maxElementId) + doRead = false; + + if (doRead) + { + if (es < bcMeta.minElementId && ee < bcMeta.minElementId) // clip lower + range_max = ee; + + if (es < bcMeta.maxElementId && ee > bcMeta.maxElementId) // clip higher + range_min = es; + + range_min = std::max(es, bcMeta.minElementId); + range_max = std::min(ee, bcMeta.maxElementId); + const cgsize_t range_num = range_max - range_min + 1; + + vertexIDs.resize(range_num * verticesPerElement, + -1234567); + + cg_elements_partial_read(cgid, base, zone, section, range_min, + range_max, vertexIDs.data(), nullptr); + + cgsize_t counter = range_min; + + for (size_t i = 0; i < vertexIDs.size(); i += verticesPerElement) + { + if (elementsToConsider.count(counter)) + { + auto last = std::min(vertexIDs.size(), i + verticesPerElement); + std::vector vec(vertexIDs.begin() + i, vertexIDs.begin() + last); + // std::cout << counter << " "; + // for (const auto &j : vec) + // std::cout << j << " "; + // std::cout << std::endl; + + for (const auto &k : vec) + { + const auto zb = k - 1; // make zero-based + auto iter = globalToVert.find(zb); + if (iter != globalToVert.end()) + bcInfo.vertexIds.insert(zb); + } + } + counter++; + } + } + } + } + else + { + vertexIDs = bcMeta.bcElementIds; + // std::cout << "Vertex BCS: "; + // for (const auto &j : vertexIDs) + // std::cout << j << " "; + // std::cout << std::endl; + for (const auto &k : vertexIDs) + { + const auto zb = k - 1; // make zero-based + auto iter = globalToVert.find(zb); + if (iter != globalToVert.end()) + bcInfo.vertexIds.insert(zb); + } + } + } +} + +apf::Mesh2 *DoIt(gmi_model *g, const std::string &fname, apf::CGNSBCMap &cgnsBCMap, const std::vector> &readMeshData) +{ + static_assert(std::is_same::value, "cgsize_t not compiled as int"); + + int cgid = -1; + auto comm = PCU_Get_Comm(); + cgp_mpi_comm(comm); + cgp_pio_mode(CGNS_ENUMV(CGP_INDEPENDENT)); + cgp_open(fname.c_str(), CGNS_ENUMV(CG_MODE_READ), &cgid); + + int nbases = -1; + cg_nbases(cgid, &nbases); + if (nbases > 1) + { + std::cout << "CGNS file has " << nbases << " nbases, only used base=1" << std::endl; + } + + std::string basename; + basename.resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + + int cellDim = -1; + int physDim = -1; + const int base = 1; + cg_base_read(cgid, base, &basename[0], &cellDim, &physDim); + const int readDim = cellDim; + + // Salome cgns is a bit on the odd side: cellDim, physDim, ncoords are not always consistent + apf::Mesh2 *mesh = apf::makeEmptyMdsMesh(g, cellDim, false); + apf::GlobalToVert globalToVert; + + LocalElementRanges localElementRanges; + using NewElements = std::vector; + std::vector localElements; + + int nzones = -1; + cg_nzones(cgid, base, &nzones); + basename = std::string(basename.c_str()); + + int nfam = -1; + cg_nfamilies(cgid, base, &nfam); + + std::vector bcInfos; + + for (int zone = 1; zone <= nzones; ++zone) + { + std::string zoneName; + zoneName.resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + + /* + 3D unstructured: NVertex, NCell3D, NBoundVertex + 2D unstructured: NVertex, NCell2D, NBoundVertex + */ + std::array sizes; + cg_zone_read(cgid, base, zone, &zoneName[0], sizes.data()); + + zoneName = std::string(zoneName.c_str()); + + int ngrids = -1; + cg_ngrids(cgid, base, zone, &ngrids); + if (ngrids > 1) + { + std::cout << __LINE__ << " CGNS is dead " << std::endl; + Kill(cgid); + } + int ncoords = -1; + cg_ncoords(cgid, base, zone, &ncoords); + + CGNS_ENUMT(ZoneType_t) + zonetype = CGNS_ENUMV(ZoneTypeNull); + cg_zone_type(cgid, base, zone, &zonetype); + + int nsections = -1; + if (zonetype == CGNS_ENUMV(Unstructured)) + { + cg_nsections(cgid, base, zone, &nsections); + } + else + { + std::cout << __LINE__ << " CGNS is dead " << std::endl; + Kill(cgid); + } + + int nBocos = -1; + cg_nbocos(cgid, base, zone, &nBocos); + + for (int section = 1; section <= nsections; section++) + { + std::string sectionName; + sectionName.resize(CGIO_MAX_NAME_LENGTH + 1, ' '); + + CGNS_ENUMT(ElementType_t) + elementType = CGNS_ENUMV(ElementTypeNull); + cgsize_t el_start = -1; + cgsize_t el_end = -1; + int num_bndry = -1; + int parent_flag = -1; + cgsize_t numElements = -1; + int verticesPerElement = -1; + + cg_section_read(cgid, base, zone, section, §ionName[0], + &elementType, &el_start, &el_end, + &num_bndry, &parent_flag); + + numElements = el_end - el_start + 1; + + cg_npe(elementType, &verticesPerElement); + + const auto readElementsAndVerts = [&](const apf::Mesh2::Type &type) { + const auto &ret = ReadElements(cgid, base, zone, section, el_start, el_end, numElements, verticesPerElement, localElementRanges); + if (std::get<1>(ret) > 0) + { + const std::vector vertexIDs = std::get<0>(ret); + localElements.emplace_back(apf::assemble(mesh, vertexIDs.data(), std::get<1>(ret), type, globalToVert)); // corresponding finalize below + const auto nverts = sizes[0]; + const auto ordinates = ReadCGNSCoords(cgid, base, zone, ncoords, nverts, vertexIDs, globalToVert); + + for (const auto &p : globalToVert) + { + const auto pp = ordinates.at(p.first); + apf::Vector3 point(pp[0], pp[1], pp[2]); + mesh->setPoint(p.second, 0, point); + + // Don't know why I wrote it like this... + // auto iter = globalToVert.find(p.first); + // if (iter != globalToVert.end()) + // { + // mesh->setPoint(iter->second, 0, point); + // } + // else + // { + // Kill(cgid, "GlobalToVert lookup problem"); + // } + } + } + }; + + if (elementType == CGNS_ENUMV(BAR_2)) + { + if (readDim == 1) + readElementsAndVerts(apf::Mesh2::EDGE); + } + else if (elementType == CGNS_ENUMV(QUAD_4)) + { + if (readDim == 2) + readElementsAndVerts(apf::Mesh2::QUAD); + } + else if (elementType == CGNS_ENUMV(TRI_3)) + { + if (readDim == 2) + readElementsAndVerts(apf::Mesh2::TRIANGLE); + } + else if (elementType == CGNS_ENUMV(TETRA_4)) + { + if (readDim == 3) + readElementsAndVerts(apf::Mesh2::TET); + } + else if (elementType == CGNS_ENUMV(PYRA_5)) + { + if (readDim == 3) + readElementsAndVerts(apf::Mesh2::PYRAMID); + } + else if (elementType == CGNS_ENUMV(HEXA_8)) + { + if (readDim == 3) + readElementsAndVerts(apf::Mesh2::HEX); + } + else + { + std::cout << __LINE__ << " CGNS is dead " + << " " << SupportedCGNSElementTypeToString(elementType) << std::endl; + Kill(cgid); + } + } + + if (nBocos > 0) + { + std::cout << std::endl; + std::cout << "Attempting to read BCS info " + << " " << nBocos << std::endl; + ReadBCInfo(cgid, base, zone, nBocos, physDim, cellDim, nsections, bcInfos, globalToVert); + std::cout << std::endl; + } + + int nsols = -1; + if (cg_nsols(cgid, base, zone, &nsols)) + Kill(cgid, "1, ", nsols); + + std::vector meshData; + if (nsols > 0) + { + for (int ns = 1; ns <= nsols; ns++) + { + CGNS_ENUMT(GridLocation_t) + location; + char sname[33]; + if (cg_sol_info(cgid, base, zone, ns, + sname, &location)) + Kill(cgid, "2"); + + int nflds = -1; + if (cg_nfields(cgid, base, zone, ns, &nflds)) + Kill(cgid, "3"); + + if (nflds > 0) + { + for (int f = 1; f <= nflds; f++) + { + CGNS_ENUMT(DataType_t) + datatype; + char name[33]; + if (cg_field_info(cgid, base, zone, ns, f, + &datatype, name)) + Kill(cgid, "4"); + + //std::cout << sname << " " << name << " " << f << " " << ns << " " << cg_DataTypeName(datatype) << " " << cg_GridLocationName(location) << std::endl; + meshData.push_back(MeshData(ns, f, datatype, location, name)); + } + } + } + } + + const auto index = [](const std::string &name) { + for (int i = 0; i < 9; i++) // check for vectors and matrices + { + // pattern matching component writer in apfCGNS.cc, patterns must be kept in-sync with this file + const std::string end("_[" + std::to_string(i) + "]"); + const std::size_t found = name.find(end); + if (found != std::string::npos) // does name contain a [x] where x[0,9) + { + return i; + } + } + return -1; + }; + + // process meshData and find vectors (and matrices) + const auto findMatch = [&meshData, &index, &cgid](MeshData &other, std::vector &meshDataGroups) { + MeshDataGroup group; + for (auto &md : meshData) + { + if (md != other) // don't compare with yourself + { + if (md.process) // make sure md has not be marked as a component already + { + if (md.name.size() == other.name.size()) + { + for (int i = 0; i < 9; i++) // check for vectors and matrices + { + // pattern matching component writer in apfCGNS.cc, patterns must be kept in-sync with this file + const std::string end("_[" + std::to_string(i) + "]"); + const std::size_t found = md.name.find(end); + if (found != std::string::npos) // does md contain a [x] where x[0,9) + { + const auto stub = md.name.substr(0, found); // get the part of the string not with a component index + const auto stubOther = other.name.substr(0, found); // get the part of the string not with a component index + + if (stub == stubOther) + { + if (md.location == other.location) // are at the same location in the mesh + { + if (md.datatype == other.datatype) // have same data type + { + md.process = false; + other.process = false; + const auto otherIndex = index(other.name); + if (otherIndex != -1) + group.insert(otherIndex, other); + else + Kill(cgid, "Bad fail"); + group.insert(i, md); + } + } + } + } + } + } + } + } + } + + if (group.size() > 0) + meshDataGroups.push_back(group); + }; + + // If a field or tag is already present, don't over-write or try to add + // Can modify this later for different behaviour, but useful default. + std::vector existingNames; + { + apf::DynamicArray tags; + mesh->getTags(tags); + for (std::size_t i = 0; i < tags.getSize(); ++i) + { + apf::MeshTag *t = tags[i]; + const std::string tagName(mesh->getTagName(t)); + existingNames.push_back(tagName); + } + for (int i = 0; i < mesh->countFields(); ++i) + { + apf::Field *f = mesh->getField(i); + const std::string fieldName(f->getName()); + existingNames.push_back(fieldName); + } + } + + for (auto &md : meshData) + { + if (md.process) + { + for (const auto &n : existingNames) + { + if (n == md.name) + md.process = false; + } + } + } + + std::vector meshDataGroups; + for (auto &md : meshData) + { + if (md.process) + findMatch(md, meshDataGroups); + } + + for (auto &md : meshData) + { + if (md.process) + { + MeshDataGroup group; + group.insert(0, md); + meshDataGroups.push_back(group); + md.process = false; + } + } + + for (auto &md : meshDataGroups) + { + if (md.process) + { + bool read = false; + for (const auto &n : readMeshData) + { + if (n.second == md.name()) + { + if (n.first == cg_GridLocationName(md.location())) + { + std::cout << md.name() << " " << n.first << " " << n.second << " " << cg_GridLocationName(md.location()) << std::endl; + read = true; + } + } + } + if (!read) + md.process = false; + } + } + + // add all double's as fields, even though they once may have been tags, since + // that info is lost in the cgns-writer (cos of cgns...) + for (const auto &md : meshDataGroups) + { + if (md.process) + { + if (md.location() == CGNS_ENUMV(Vertex)) + { + if (md.datatype() == CGNS_ENUMV(Integer)) + { + } + else if (md.datatype() == CGNS_ENUMV(RealDouble)) + { + const cgsize_t lowest = globalToVert.begin()->first + 1; // one based + const cgsize_t highest = globalToVert.rbegin()->first + 1; // one based + + cgsize_t range_min[3]; + range_min[0] = lowest; + range_min[1] = lowest; + range_min[2] = lowest; + cgsize_t range_max[3]; + range_max[0] = highest; + range_max[1] = highest; + range_max[2] = highest; + const cgsize_t numToRead = range_max[0] - range_min[0] + 1; // one based + + apf::Field *field = nullptr; + if (md.size() == 1) + field = apf::createFieldOn(mesh, md.name().c_str(), apf::SCALAR); + else if (md.size() == 3) + field = apf::createFieldOn(mesh, md.name().c_str(), apf::VECTOR); + else if (md.size() == 9) + field = apf::createFieldOn(mesh, md.name().c_str(), apf::MATRIX); + else + Kill(cgid, "Tensor size not accounted for"); + + double scalar(-123456.0); + apf::Vector3 vector3(-123456.0, -123456.0, -123456.0); + apf::Matrix3x3 matrix3x3(-123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0); + + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = mesh->begin(0); + while ((elem = mesh->iterate(it))) + { + if (md.size() == 1) + apf::setScalar(field, elem, 0, scalar); + else if (md.size() == 3) + apf::setVector(field, elem, 0, vector3); + else if (md.size() == 9) + apf::setMatrix(field, elem, 0, matrix3x3); + else + Kill(cgid, "Tensor size not accounted for"); + } + mesh->end(it); + + using CGNSType = double; + std::vector meshVals(numToRead, -123456.0); + for (std::size_t i = 0; i < md.size(); i++) + { + if (cgp_field_read_data(cgid, base, zone, md.sIndex(i), md.fIndex(i), + &range_min[0], &range_max[0], meshVals.data())) + Kill(cgid, "Failed cgp_field_read_data"); + + cgsize_t counter = lowest; + for (cgsize_t it = 0; it < numToRead; it++) + { + cgsize_t zeroBased = counter - 1; // remove as per the addition above + auto iter = globalToVert.find(zeroBased); + if (iter != globalToVert.end()) + { + auto *elem = iter->second; + if (md.size() == 1) + { + apf::setScalar(field, elem, 0, meshVals.at(it)); + } + else if (md.size() == 3) + { + apf::getVector(field, elem, 0, vector3); + vector3[i] = meshVals.at(it); + apf::setVector(field, elem, 0, vector3); + } + else if (md.size() == 9) + { + apf::getMatrix(field, elem, 0, matrix3x3); + if (i < 3) + matrix3x3[0][i] = meshVals.at(it); + else if (i > 2 && i < 6) + matrix3x3[1][i - 3] = meshVals.at(it); + else if (i > 5 && i < 9) + matrix3x3[2][i - 6] = meshVals.at(it); + apf::setMatrix(field, elem, 0, matrix3x3); + } + else + Kill(cgid, "Tensor size not accounted for"); + } + counter++; + } + } + } + else if (md.datatype() == CGNS_ENUMV(LongInteger)) + { + } + else + { + Kill(cgid, "Don't know how to process this at the moment"); + } + } + else if (md.location() == CGNS_ENUMV(CellCenter)) + { + if (md.datatype() == CGNS_ENUMV(Integer)) + { + } + else if (md.datatype() == CGNS_ENUMV(RealDouble)) + { + // HUGE assumption here, I assume the order the cgns elements are read (and therefore their global number in the file) + // is the same order they are created and added to the mesh by buildElements. + // I combine localElementRanges to give me the global indices, with vector and hope it works + PCU_ALWAYS_ASSERT_VERBOSE(localElementRanges.size() == localElements.size(), + "Size don't match for element/number ranges"); + + const auto dim = mesh->getDimension(); + apf::Field *field = nullptr; + md.info(); + if (md.size() == 1) + field = apf::createField(mesh, md.name().c_str(), apf::SCALAR, apf::getConstant(dim)); + else if (md.size() == 3) + field = apf::createField(mesh, md.name().c_str(), apf::VECTOR, apf::getConstant(dim)); + else if (md.size() == 9) + field = apf::createField(mesh, md.name().c_str(), apf::MATRIX, apf::getConstant(dim)); + else + Kill(cgid, "Tensor size not accounted for"); + + double scalar(-123456.0); + apf::Vector3 vector3(-123456.0, -123456.0, -123456.0); + apf::Matrix3x3 matrix3x3(-123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0, -123456.0); + + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = mesh->begin(dim); + while ((elem = mesh->iterate(it))) + { + if (md.size() == 1) + apf::setScalar(field, elem, 0, scalar); + else if (md.size() == 3) + apf::setVector(field, elem, 0, vector3); + else if (md.size() == 9) + apf::setMatrix(field, elem, 0, matrix3x3); + else + Kill(cgid, "Tensor size not accounted for"); + } + mesh->end(it); + + for (std::size_t r = 0; r < localElementRanges.size(); r++) + { + const auto range = localElementRanges[r]; + const cgsize_t lowest = range.first; // one based + const cgsize_t highest = range.second; // one based + + cgsize_t range_min[3]; + range_min[0] = lowest; + range_min[1] = lowest; + range_min[2] = lowest; + cgsize_t range_max[3]; + range_max[0] = highest; + range_max[1] = highest; + range_max[2] = highest; + const std::size_t numToRead = range_max[0] - range_min[0] + 1; // one based + PCU_ALWAYS_ASSERT_VERBOSE(numToRead == localElements[r].size(), + "Size don't match for element/number sub-ranges"); + + using CGNSType = double; + std::vector meshVals(numToRead, -123456.0); + for (std::size_t i = 0; i < md.size(); i++) + { + if (cgp_field_read_data(cgid, base, zone, md.sIndex(i), md.fIndex(i), + &range_min[0], &range_max[0], meshVals.data())) + Kill(cgid, "Failed cgp_field_read_data"); + + for (std::size_t it = 0; it < numToRead; it++) + { + elem = localElements[r][it]; + + if (md.size() == 1) + { + apf::setScalar(field, elem, 0, meshVals.at(it)); + } + else if (md.size() == 3) + { + apf::getVector(field, elem, 0, vector3); + vector3[i] = meshVals.at(it); + apf::setVector(field, elem, 0, vector3); + } + else if (md.size() == 9) + { + apf::getMatrix(field, elem, 0, matrix3x3); + if (i < 3) + matrix3x3[0][i] = meshVals.at(it); + else if (i > 2 && i < 6) + matrix3x3[1][i - 3] = meshVals.at(it); + else if (i > 5 && i < 9) + matrix3x3[2][i - 6] = meshVals.at(it); + apf::setMatrix(field, elem, 0, matrix3x3); + } + else + { + Kill(cgid, "Tensor size not accounted for"); + } + } + } + } + } + else if (md.datatype() == CGNS_ENUMV(LongInteger)) + { + } + else + { + Kill(cgid, "Don't know how to process this at the moment"); + } + } + else + { + Kill(cgid, "Don't know how to process this at the moment"); + } + } + } + } + + // free up memory + if (PCU_Comm_Initialized()) + cgp_close(cgid); + else + cg_close(cgid); + + { + apf::MeshTag *tag = nullptr; + tag = mesh->createIntTag("origCGNSGlobalVertID", 1); // 1 is size of tag + const cgsize_t lowest = globalToVert.begin()->first + 1; // one based + const cgsize_t highest = globalToVert.rbegin()->first + 1; // one based + const cgsize_t numToRead = highest - lowest + 1; // one based + cgsize_t counter = lowest; + for (cgsize_t it = 0; it < numToRead; it++) + { + cgsize_t zeroBased = counter - 1; // remove as per the addition above + auto iter = globalToVert.find(zeroBased); + if (iter != globalToVert.end()) + { + const int oneBased = zeroBased + 1; + mesh->setIntTag(iter->second, tag, &oneBased); + } + counter++; + } + } + + { + apf::MeshTag *tag = nullptr; + tag = mesh->createIntTag("origCGNSGlobalElemID", 1); // 1 is size of tag + for (std::size_t r = 0; r < localElementRanges.size(); r++) + { + const auto range = localElementRanges[r]; + const cgsize_t lowest = range.first; // one based + const cgsize_t highest = range.second; // one based + const std::size_t numToRead = highest - lowest + 1; // one based + int counter = lowest; + for (std::size_t it = 0; it < numToRead; it++) + { + apf::MeshEntity *elem = nullptr; + elem = localElements[r][it]; + mesh->setIntTag(elem, tag, &counter); + counter++; + } + } + } + // not sure of the order of these three + // and with reference to: + + apf::finalise(mesh, globalToVert); + mesh->acceptChanges(); + apf::alignMdsRemotes(mesh); + + { + apf::GlobalNumbering *gn = nullptr; + gn = apf::makeGlobal(apf::numberOwnedNodes(mesh, "vert Idx")); + apf::synchronize(gn); + } + // no synchronize call + // https://github.com/SNLComputation/Albany/blob/master/src/disc/pumi/Albany_APFDiscretization.cpp @ various place throughout file + // https://github.com/SCOREC/core/issues/249 + { + apf::makeGlobal(apf::numberElements(mesh, "elem Idx")); + } + + for (auto &bc : bcInfos) + { + bc.TagVertices(cgid, mesh, globalToVert); + } + + for (auto &bc : bcInfos) + { + bc.TagBCEntities(cgid, mesh, cgnsBCMap); + } + + apf::deriveMdsModel(mesh); + apf::verify(mesh, true); + + return mesh; +} // namespace + +apf::Mesh2 *DoIt(gmi_model *g, const std::string &fname, apf::CGNSBCMap &cgnsBCMap) +{ + std::vector> meshData; + return DoIt(g, fname, cgnsBCMap, meshData); +} + +} // namespace + +namespace apf +{ + +// caller needs to bring up and pull down mpi/pcu: mpi/pcu is required and assumed. +Mesh2 *loadMdsFromCGNS(gmi_model *g, const char *fname, apf::CGNSBCMap &cgnsBCMap, const std::vector> &meshData) +{ +#ifdef HAVE_CGNS + Mesh2 *m = DoIt(g, fname, cgnsBCMap, meshData); + return m; +#else + Mesh2 *m = nullptr; + PCU_ALWAYS_ASSERT_VERBOSE(m != nullptr, + "Build with ENABLE_CGNS to allow this functionality."); + exit(EXIT_FAILURE); + return m; +#endif +} + +// caller needs to bring up and pull down mpi/pcu: mpi/pcu is required and assumed. +Mesh2 *loadMdsFromCGNS(gmi_model *g, const char *fname, apf::CGNSBCMap &cgnsBCMap) +{ +#ifdef HAVE_CGNS + Mesh2 *m = DoIt(g, fname, cgnsBCMap); + return m; +#else + Mesh2 *m = nullptr; + PCU_ALWAYS_ASSERT_VERBOSE(m != nullptr, + "Build with ENABLE_CGNS to allow this functionality."); + exit(EXIT_FAILURE); + return m; +#endif +} + +} // namespace apf diff --git a/mds/pkg_tribits.cmake b/mds/pkg_tribits.cmake index b885dc256..07211c87b 100644 --- a/mds/pkg_tribits.cmake +++ b/mds/pkg_tribits.cmake @@ -22,6 +22,10 @@ set(MDS_SOURCES mdsGmsh.cc mdsUgrid.cc) +if(ENABLE_CGNS) + set(MDS_SOURCES ${MDS_SOURCES} mdsCGNS.cc) +endif(ENABLE_CGNS) + set(MDS_HEADERS apfMDS.h apfBox.h diff --git a/omega_h/CMakeLists.txt b/omega_h/CMakeLists.txt index 9a80399a1..750fa0538 100644 --- a/omega_h/CMakeLists.txt +++ b/omega_h/CMakeLists.txt @@ -17,7 +17,11 @@ set(HEADERS # Add the apf_omega_h library add_library(apf_omega_h ${SOURCES}) -target_compile_options(apf_omega_h PUBLIC "-std=c++11") +if(ENABLE_CGNS) + target_compile_options(apf_omega_h PUBLIC "-std=c++14") +else() + target_compile_options(apf_omega_h PUBLIC "-std=c++11") +endif() # Include directories target_include_directories(apf_omega_h INTERFACE diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a19e4fd56..709492f91 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -43,6 +43,9 @@ util_exe_func(extrude extrude.cc) # Mesh conversion/formatting utilities util_exe_func(from_gmsh gmsh.cc) +if(ENABLE_CGNS) + util_exe_func(from_cgns cgns.cc) +endif() util_exe_func(from_ansys ansys.cc) util_exe_func(from_neper neper.cc) util_exe_func(from_ugrid ugrid.cc) diff --git a/test/cgns.cc b/test/cgns.cc new file mode 100644 index 000000000..6f853d2e0 --- /dev/null +++ b/test/cgns.cc @@ -0,0 +1,553 @@ +/* Author: Dr Andrew Parker (2019) - FGE Ltd */ + +#include +#include +#include +#include +#include +#include +#include +#include +// +#include +#include +// +#include +#include +#include +#include + +// https://gist.github.com/bgranzow/98087114166956646da684ed98acab02 +apf::MeshTag *create_int_tag(const std::string &name, apf::Mesh2 *m, int dim) +{ + apf::MeshTag *tag = m->createIntTag(name.c_str(), 1); // 1 is size of tag + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(dim); + int vals[1]; + vals[0] = PCU_Comm_Self(); + while ((elem = m->iterate(it))) + m->setIntTag(elem, tag, vals); + m->end(it); + return tag; +} + +//https://github.com/CEED/PUMI/blob/master/ma/maDBG.cc +apf::Field *convert_tag_doubleField(const std::string &name, apf::Mesh2 *m, apf::MeshTag *t, int dim) +{ + apf::MeshEntity *elem = nullptr; + apf::MeshIterator *it = m->begin(dim); + apf::Field *f = nullptr; + const auto fieldName = name.c_str(); + if (dim == 0) + f = apf::createFieldOn(m, fieldName, apf::SCALAR); + else + f = apf::createField(m, fieldName, apf::SCALAR, apf::getConstant(dim)); + + int vals[1]; + while ((elem = m->iterate(it))) + { + m->getIntTag(elem, t, vals); + double dval; + dval = vals[0]; + apf::setScalar(f, elem, 0, dval); + } + m->end(it); + return f; +} + +void balance(const std::string &prefix, const apf::ZoltanMethod &method, apf::Mesh2 *m) +{ + const auto dim = m->getDimension(); + convert_tag_doubleField("procID_prebalance", m, create_int_tag("procID_prebalance", m, dim), dim); + + apf::MeshTag *weights = Parma_WeighByMemory(m); + apf::Balancer *balancer = makeZoltanBalancer(m, method, apf::REPARTITION); + balancer->balance(weights, 1.10); + delete balancer; + apf::removeTagFromDimension(m, weights, m->getDimension()); + Parma_PrintPtnStats(m, ""); + m->destroyTag(weights); + + const std::string name = prefix + "_balance_" + std::to_string(PCU_Comm_Peers()) + "procs"; + apf::writeVtkFiles(name.c_str(), m); +} + +//https://github.com/SCOREC/core/blob/4b854ae996cb261a22f6ee6b704569b78866004c/test/reorder.cc +void simpleReorder(const std::string &prefix, apf::Mesh2 *m) +{ + { + apf::GlobalNumbering *gn = nullptr; + gn = apf::makeGlobal(apf::numberOwnedNodes(m, "vert Idx_preOrder")); + apf::synchronize(gn); + } + // no synchronize call + // https://github.com/SNLComputation/Albany/blob/master/src/disc/pumi/Albany_APFDiscretization.cpp @ various place throughout file + // https://github.com/SCOREC/core/issues/249 + { + apf::makeGlobal(apf::numberElements(m, "elem Idx_preOrder")); + } + + //apf::MeshTag *order = Parma_BfsReorder(m); + //apf::reorderMdsMesh(m, order); + apf::reorderMdsMesh(m); + + { + apf::GlobalNumbering *gn = nullptr; + gn = apf::makeGlobal(apf::numberOwnedNodes(m, "vert Idx_pstOrder")); + apf::synchronize(gn); + } + // no synchronize call + // https://github.com/SNLComputation/Albany/blob/master/src/disc/pumi/Albany_APFDiscretization.cpp @ various place throughout file + // https://github.com/SCOREC/core/issues/249 + { + apf::makeGlobal(apf::numberElements(m, "elem Idx_pstOrder")); + } + + const std::string name = prefix + "_reorder_" + std::to_string(PCU_Comm_Peers()) + "procs"; + apf::writeVtkFiles(name.c_str(), m); +} + +pMesh toPumi(const std::string &prefix, gmi_model *g, apf::Mesh2 *mesh) +{ + //create the pumi instance + pumi::instance()->model = new gModel(g); + pMesh pm = pumi_mesh_load(mesh); + std::cout << pm << std::endl; + pumi_mesh_verify(pm); + const std::string name = prefix + "_toPUMI_" + std::to_string(PCU_Comm_Peers()) + "procs"; + pumi_mesh_write(pm, name.c_str(), "vtk"); + return pm; +} + +auto additional(const std::string &prefix, gmi_model *g, apf::Mesh2 *mesh) +{ + // seems essential to make pm first before calling balance or reorder... + auto pm = toPumi(prefix, g, mesh); + balance(prefix, apf::RCB, pm); + simpleReorder(prefix, pm); + + //create an element field + const int mdim = pumi_mesh_getDim(pm); + pShape s = pumi_shape_getConstant(mdim); + const int dofPerElm = 1; + pField f = pumi_field_create(pm, "elmField", dofPerElm, PUMI_PACKED, s); + + pMeshIter it = pm->begin(mdim); + pMeshEnt e; + double v = 0; + while ((e = pm->iterate(it))) + pumi_node_setField(f, e, 0, &v); + pm->end(it); + + const int ghost = mdim; + const int bridge = ghost - 1; + const int numLayers = 1; + const int ghostAcrossCopies = 1; + pumi_ghost_createLayer(pm, bridge, ghost, numLayers, ghostAcrossCopies); + + it = pm->begin(mdim); + v = 1; + while ((e = pm->iterate(it))) + { + if (!pumi_ment_isGhost(e)) + pumi_node_setField(f, e, 0, &v); + } + pm->end(it); + + const std::string name = prefix + "_toPUMI_" + std::to_string(PCU_Comm_Peers()) + "procs"; + // only the owned elements will have a elmField value of 1 + pumi_mesh_write(pm, (name + "_beforeSync").c_str(), "vtk"); + + pumi_field_synchronize(f); + + // owned and ghosted elements will have a elmField value of 1 + pumi_mesh_write(pm, (name + "_afterSync").c_str(), "vtk"); + + const auto clean = [&pm, &f]() { + // clean-up + pumi_field_delete(f); + pumi_ghost_delete(pm); + pumi_mesh_delete(pm); + }; + return clean; +} + +std::string doit(apf::CGNSBCMap &cgnsBCMap, const std::string &argv1, const std::string &argv2, const std::string &post, const bool &additionalTests, const bool &writeCGNS, const std::vector> &meshData) +{ + gmi_register_null(); + gmi_register_mesh(); + gmi_model *g = gmi_load(".null"); + apf::Mesh2 *m = apf::loadMdsFromCGNS(g, argv1.c_str(), cgnsBCMap, meshData); + m->verify(); + // + m->writeNative(argv2.c_str()); + // so we can see the result + const std::string path = argv1.c_str(); + const std::size_t found = path.find_last_of("/\\"); + std::size_t index = found + 1; + if (found == std::string::npos) + index = 0; + + //std::cout << "FOUND " << found << " " << path.size() << " " << path << " " << index << std::endl; + +#ifndef NDEBUG // debug settings, cmake double negative.... + const auto prefix = path.substr(index) + "_debug" + post; +#else // optimised setting + const auto prefix = path.substr(index) + "_release" + post; +#endif + + const auto dim = m->getDimension(); + if (dim == 3) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh"); + apf::writeVtkFiles(name.c_str(), m, 3); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_faceMesh"); + apf::writeVtkFiles(name.c_str(), m, 2); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_edgeMesh"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + else if (dim == 2) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh"); + apf::writeVtkFiles(name.c_str(), m, 2); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_edgeMesh"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + else if (dim == 1) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + + std::string cgnsOutputName; + if (writeCGNS) + { + cgnsOutputName = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + "_outputFile.cgns"; + apf::writeCGNS(cgnsOutputName.c_str(), m, cgnsBCMap); + } + + // main purpose is to call additional tests through the test harness testing.cmake + std::function cleanUp; + if (additionalTests) + { + cleanUp = additional(prefix, g, m); + // + } + + if (additionalTests) + { + // add dummy matrix to mesh + const auto addMatrix = [](apf::Mesh2 *mesh, const int &dim) { + apf::Field *field = nullptr; + const std::string name = "DummyMatrix_" + std::to_string(dim); + apf::MeshTag *tag = nullptr; + if (dim != 0) + { + field = apf::createField(mesh, name.c_str(), apf::MATRIX, apf::getConstant(dim)); + if (dim == mesh->getDimension()) + tag = mesh->findTag("origCGNSGlobalElemID"); + } + else if (dim == 0) + { + field = apf::createFieldOn(mesh, name.c_str(), apf::MATRIX); + tag = mesh->findTag("origCGNSGlobalVertID"); + } + //std::cout << "*************************** "<< dim << std::endl; + apf::MeshIterator *it = mesh->begin(dim); + apf::MeshEntity *elm = nullptr; + int vals[1]; + while ((elm = mesh->iterate(it))) + { + apf::Matrix3x3 m(11, 12, 13, 21, 22, 23, 31, 32, 33); + if (tag) + { + mesh->getIntTag(elm, tag, vals); + m[0][0] = vals[0]; + m[1][1] = -vals[0]; + m[2][2] = vals[0] * vals[0]; + } + apf::setMatrix(field, elm, 0, m); + } + mesh->end(it); + //apf::destroyField(field); + }; + + for (int i = 0; i <= m->getDimension(); i++) + addMatrix(m, i); + + // add dummy vector to mesh + const auto addVector = [](apf::Mesh2 *mesh, const int &dim) { + apf::Field *field = nullptr; + const std::string name = "DummyVector_" + std::to_string(dim); + apf::MeshTag *tag = nullptr; + if (dim != 0) + { + field = apf::createField(mesh, name.c_str(), apf::VECTOR, apf::getConstant(dim)); + if (dim == mesh->getDimension()) + tag = mesh->findTag("origCGNSGlobalElemID"); + } + else if (dim == 0) + { + field = apf::createFieldOn(mesh, name.c_str(), apf::VECTOR); + tag = mesh->findTag("origCGNSGlobalVertID"); + } + //std::cout << "*************************** "<< dim << std::endl; + apf::MeshIterator *it = mesh->begin(dim); + apf::MeshEntity *elm = nullptr; + int vals[1]; + while ((elm = mesh->iterate(it))) + { + apf::Vector3 v; + if (tag) + { + mesh->getIntTag(elm, tag, vals); + v[0] = vals[0]; + v[1] = -vals[0]; + v[2] = vals[0] * vals[0]; + } + else + { + v[0] = 1.0; + v[1] = 2.0; + v[2] = 3.0; + } + + apf::setVector(field, elm, 0, v); + } + mesh->end(it); + //apf::destroyField(field); + }; + + for (int i = 0; i <= m->getDimension(); i++) + addVector(m, i); + + // add dummy scalar to mesh + const auto addScalar = [](apf::Mesh2 *mesh, const int &dim) { + apf::Field *field = nullptr; + const std::string name = "DummyScalar_" + std::to_string(dim); + apf::MeshTag *tag = nullptr; + if (dim != 0) + { + field = apf::createField(mesh, name.c_str(), apf::SCALAR, apf::getConstant(dim)); + if (dim == mesh->getDimension()) + tag = mesh->findTag("origCGNSGlobalElemID"); + } + else if (dim == 0) + { + field = apf::createFieldOn(mesh, name.c_str(), apf::SCALAR); + tag = mesh->findTag("origCGNSGlobalVertID"); + } + //std::cout << "*************************** "<< dim << std::endl; + apf::MeshIterator *it = mesh->begin(dim); + apf::MeshEntity *elm = nullptr; + int vals[1]; + vals[0] = -1234567; + while ((elm = mesh->iterate(it))) + { + double v = 1.0; + if (tag) + { + mesh->getIntTag(elm, tag, vals); + v = double(vals[0]); + } + apf::setScalar(field, elm, 0, v); + } + mesh->end(it); + //apf::destroyField(field); + }; + + for (int i = 0; i <= m->getDimension(); i++) + addScalar(m, i); + + if (dim == 3) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 3); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_faceMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 2); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_edgeMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + else if (dim == 2) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 2); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_edgeMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + else if (dim == 1) + { + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_cellMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 1); + } + { + const auto name = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + std::string("_toVTK_vertexMesh_additional"); + apf::writeVtkFiles(name.c_str(), m, 0); + } + } + + if (writeCGNS) + { + // what this one to be re-read if doing re-reading so that the dummy variables (vector/matrix/scalar) are there + cgnsOutputName = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + "_additional_outputFile.cgns"; + apf::writeCGNS(cgnsOutputName.c_str(), m, cgnsBCMap); + } + // + { + apf::Field *field = nullptr; + const std::string name = "displace"; + field = apf::createFieldOn(m, name.c_str(), apf::VECTOR); + apf::MeshIterator *it = m->begin(0); + apf::MeshEntity *elm = nullptr; + while ((elm = m->iterate(it))) + { + apf::Vector3 v; + v[0] = 0.25; + v[1] = 0.25; + v[2] = 0.35; + apf::setVector(field, elm, 0, v); + } + m->end(it); + apf::displaceMesh(m, field); + if (writeCGNS) + { + std::string cgnsOutputName = prefix + "_" + std::to_string(PCU_Comm_Peers()) + "procs" + "_additional_moved_outputFile.cgns"; + apf::writeCGNS(cgnsOutputName.c_str(), m, cgnsBCMap); + } + } + } + + if (additionalTests) + { + //cleanUp(); + } + else if (!additionalTests) + { + m->destroyNative(); + apf::destroyMesh(m); + } + // + return cgnsOutputName; +} + +int main(int argc, char **argv) +{ +#ifdef HAVE_CGNS + MPI_Init(&argc, &argv); + PCU_Comm_Init(); + lion_set_verbosity(1); + bool additionalTests = false; + if (argc < 3) + { + if (!PCU_Comm_Self()) + printf("Usage: %s \n", argv[0]); + MPI_Finalize(); + exit(EXIT_FAILURE); + return -1; + } + else if (argc == 4) + { + if (std::string(argv[3]) == "additional") + additionalTests = true; + else + { + if (!PCU_Comm_Self()) + printf("Usage: %s additional\n", argv[0]); + MPI_Finalize(); + exit(EXIT_FAILURE); + return -1; + } + } + else if (argc > 4) + { + if (!PCU_Comm_Self()) + printf("Usage: %s \n", argv[0]); + MPI_Finalize(); + exit(EXIT_FAILURE); + return -1; + } + + // Phase 1 + std::string cgnsOutputName; + { + apf::CGNSBCMap cgnsBCMap; + std::vector> meshData; + cgnsOutputName = doit(cgnsBCMap, argv[1], argv[2], "", additionalTests, additionalTests, meshData); + } + // Phase 2 + if (additionalTests) + { + std::vector> meshData; + meshData.push_back(std::make_pair("Vertex", "DummyScalar_0")); + meshData.push_back(std::make_pair("Vertex", "DummyVector_0")); + meshData.push_back(std::make_pair("Vertex", "DummyMatrix_0")); + + meshData.push_back(std::make_pair("CellCenter", "DummyScalar_3")); + meshData.push_back(std::make_pair("CellCenter", "DummyVector_3")); + meshData.push_back(std::make_pair("CellCenter", "DummyMatrix_3")); + + meshData.push_back(std::make_pair("CellCenter", "DummyScalar_2")); + meshData.push_back(std::make_pair("CellCenter", "DummyVector_2")); + meshData.push_back(std::make_pair("CellCenter", "DummyMatrix_2")); + + meshData.push_back(std::make_pair("CellCenter", "DummyScalar_1")); + meshData.push_back(std::make_pair("CellCenter", "DummyVector_1")); + meshData.push_back(std::make_pair("CellCenter", "DummyMatrix_1")); + apf::CGNSBCMap cgnsBCMap; + std::cout << "RE-READING " << std::endl; + doit(cgnsBCMap, cgnsOutputName.c_str(), "tempy.smb", "_reread", false, true, meshData); + } + // + PCU_Comm_Free(); + MPI_Finalize(); + return 0; +#else + PCU_ALWAYS_ASSERT_VERBOSE(true == false, + "Build with ENABLE_CGNS to allow this functionality."); + exit(EXIT_FAILURE); + return -1; +#endif +} diff --git a/test/constructThenGhost.cc b/test/constructThenGhost.cc index 4caafe06a..56cdb85e2 100644 --- a/test/constructThenGhost.cc +++ b/test/constructThenGhost.cc @@ -59,7 +59,7 @@ int main(int argc, char** argv) pMeshIter it = pm->begin(mdim); pMeshEnt e; double v = 0; - while ((e = m->iterate(it))) + while ((e = pm->iterate(it))) pumi_node_setField(f,e,0,&v); pm->end(it); @@ -71,7 +71,7 @@ int main(int argc, char** argv) it = pm->begin(mdim); v = 1; - while ((e = m->iterate(it))) { + while ((e = pm->iterate(it))) { if (!pumi_ment_isGhost(e)) pumi_node_setField(f,e,0,&v); } diff --git a/test/testing.cmake b/test/testing.cmake index e0d368fee..96c419933 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -450,6 +450,94 @@ if(ENABLE_ZOLTAN) "4" "rib" "reptn" "1" ) endif() + +if(ENABLE_CGNS) +# +# sort of an arbitrary choice +set(numProcs 4) +# +set(CGNSDIR ${MESHES}/cgns/basic) +# +# 2D tests including for mixed cells +# +mpi_test(cgns_2d_1 ${numProcs} + ./from_cgns + "${CGNSDIR}/2D/4quads.cgns" + 4quads.smb + additional) +mpi_test(cgns_2d_2 ${numProcs} + ./from_cgns + "${CGNSDIR}/2D/5quad1Tri.cgns" + 5quad1Tri.smb + additional) +mpi_test(cgns_2d_3 ${numProcs} + ./from_cgns + "${CGNSDIR}/2D/5quad2Tri.cgns" + 5quad2Tri.smb + additional) +mpi_test(cgns_2d_4 ${numProcs} + ./from_cgns + "${CGNSDIR}/2D/9tris.cgns" + 9tris.smb + additional) +# +# 3D tests including for mixed cells +# +mpi_test(cgns_3d_1 ${numProcs} + ./from_cgns + "${CGNSDIR}/3D/tets_pyra.cgns" + tets_pyra.smb + additional) +mpi_test(cgns_3d_2 ${numProcs} + ./from_cgns + "${CGNSDIR}/3D/hexs.cgns" + hexs.smb + additional) +# +# 3D BCS tests +# +set(numProcs 5) +# +set(CGNSDIR ${MESHES}/cgns/withBCS/3D) +# +mpi_test(cgns_bcs_1 ${numProcs} + ./from_cgns + "${CGNSDIR}/mixed.cgns" + bcs1.smb + additional) + +mpi_test(cgns_bcs_hex ${numProcs} + ./from_cgns + "${CGNSDIR}/8hexs.cgns" + bcshex.smb + additional) +# +# 2D BCS tests +# +set(numProcs 4) +# +set(CGNSDIR ${MESHES}/cgns/withBCS/2D) +# +mpi_test(cgns_bcs_2 ${numProcs} + ./from_cgns + "${CGNSDIR}/4quads.cgns" + bcs2.smb + additional) +# +# 1D BCS tests +# +set(numProcs 3) +# +set(CGNSDIR ${MESHES}/cgns/withBCS/1D) +# +mpi_test(cgns_bcs_3 ${numProcs} + ./from_cgns + "${CGNSDIR}/edges.cgns" + bcs3.smb + additional) + +endif(ENABLE_CGNS) + mpi_test(construct 4 ./construct "${MDIR}/cube.dmg"