From 22c1e47bf40c89f18eed2c566e2119b6b6533c18 Mon Sep 17 00:00:00 2001 From: Keisuke Inoue Date: Fri, 8 Jan 2021 13:43:11 +0900 Subject: [PATCH] Close #72 --- CMakeLists.txt | 21 ++ iriclib.pro | 21 +- iricsolver_ftoc.c | 105 +++++++++ iricsolverlib.cpp | 125 ++++++++++ iricsolverlib.h | 22 ++ iricsolverlib_api.h | 16 ++ iricsolverlib_cell2d.cpp | 55 +++++ iricsolverlib_cell2d.h | 47 ++++ iricsolverlib_grid2d.cpp | 311 +++++++++++++++++++++++++ iricsolverlib_grid2d.h | 51 ++++ iricsolverlib_point2d.h | 30 +++ iricsolverlib_quadcell.cpp | 75 ++++++ iricsolverlib_quadcell.h | 21 ++ iricsolverlib_rect2d.cpp | 120 ++++++++++ iricsolverlib_rect2d.h | 47 ++++ iricsolverlib_tricell.cpp | 64 +++++ iricsolverlib_tricell.h | 24 ++ private/iricsolverlib_cell2d_impl.h | 22 ++ private/iricsolverlib_grid2d_impl.h | 38 +++ private/iricsolverlib_point2d_detail.h | 43 ++++ 20 files changed, 1257 insertions(+), 1 deletion(-) create mode 100644 iricsolver_ftoc.c create mode 100644 iricsolverlib.cpp create mode 100644 iricsolverlib.h create mode 100644 iricsolverlib_api.h create mode 100644 iricsolverlib_cell2d.cpp create mode 100644 iricsolverlib_cell2d.h create mode 100644 iricsolverlib_grid2d.cpp create mode 100644 iricsolverlib_grid2d.h create mode 100644 iricsolverlib_point2d.h create mode 100644 iricsolverlib_quadcell.cpp create mode 100644 iricsolverlib_quadcell.h create mode 100644 iricsolverlib_rect2d.cpp create mode 100644 iricsolverlib_rect2d.h create mode 100644 iricsolverlib_tricell.cpp create mode 100644 iricsolverlib_tricell.h create mode 100644 private/iricsolverlib_cell2d_impl.h create mode 100644 private/iricsolverlib_grid2d_impl.h create mode 100644 private/iricsolverlib_point2d_detail.h diff --git a/CMakeLists.txt b/CMakeLists.txt index dec8c4b..cf986e8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,6 +40,7 @@ if(MSVC_VERSION EQUAL 1400 OR MSVC_VERSION GREATER 1400) endif() add_definitions(-DIRICLIBDLL_LIBRARY) +add_definitions(-DIRICSOLVERLIBDLL_LIBRARY) SET(iriclib_SOURCES error_macros.h @@ -71,6 +72,21 @@ iriclib_polyline.h iriclib_riversurvey.cpp iriclib_riversurvey.h iriclib_single.c +iricsolverlib.cpp +iricsolverlib.h +iricsolverlib_api.h +iricsolverlib_cell2d.cpp +iricsolverlib_cell2d.h +iricsolverlib_grid2d.cpp +iricsolverlib_grid2d.h +iricsolverlib_point2d.h +iricsolverlib_quadcell.cpp +iricsolverlib_quadcell.h +iricsolverlib_rect2d.cpp +iricsolverlib_rect2d.h +iricsolverlib_tricell.cpp +iricsolverlib_tricell.h +iricsolver_ftoc.c private/filelocker_impl.h private/iriclib_cgnsfile_baseiterativet.h private/iriclib_cgnsfile_baseiterativet_detail.h @@ -81,6 +97,9 @@ private/iriclib_cgnsfile_solutionwriterdividesolutions.cpp private/iriclib_cgnsfile_solutionwriterdividesolutions.h private/iriclib_cgnsfile_solutionwriterstandard.cpp private/iriclib_cgnsfile_solutionwriterstandard.h +private/iricsolverlib_cell2d_impl.h +private/iricsolverlib_grid2d_impl.h +private/iricsolverlib_point2d_detail.h ) # cgns include (set CMAKE_PREFIX_PATH) @@ -144,6 +163,8 @@ ${PROJECT_SOURCE_DIR}/iriclib_pointmap.h ${PROJECT_SOURCE_DIR}/iriclib_polygon.h ${PROJECT_SOURCE_DIR}/iriclib_polyline.h ${PROJECT_SOURCE_DIR}/iriclib_riversurvey.h +${PROJECT_SOURCE_DIR}/iricsolverlib.h +${PROJECT_SOURCE_DIR}/iricsolverlib_api.h ) SET(iriclib_private_Headers diff --git a/iriclib.pro b/iriclib.pro index f7dd7f1..3b3e0e0 100644 --- a/iriclib.pro +++ b/iriclib.pro @@ -14,6 +14,7 @@ unix { CONFIG += dll DEFINES += IRICLIBDLL_LIBRARY +DEFINES += IRICSOLVERLIBDLL_LIBRARY QT = DEFINES += UPPERCASE @@ -51,13 +52,24 @@ HEADERS += error_macros.h \ iriclib_polygon.h \ iriclib_polyline.h \ iriclib_riversurvey.h \ + iricsolverlib.h \ + iricsolverlib_api.h \ + iricsolverlib_cell2d.h \ + iricsolverlib_grid2d.h \ + iricsolverlib_point2d.h \ + iricsolverlib_quadcell.h \ + iricsolverlib_rect2d.h \ + iricsolverlib_tricell.h \ private/filelocker_impl.h \ private/iriclib_cgnsfile_baseiterativet.h \ private/iriclib_cgnsfile_baseiterativet_detail.h \ private/iriclib_cgnsfile_impl.h \ private/iriclib_cgnsfile_solutionwriter.h \ private/iriclib_cgnsfile_solutionwriterdividesolutions.h \ - private/iriclib_cgnsfile_solutionwriterstandard.h + private/iriclib_cgnsfile_solutionwriterstandard.h \ + private/iricsolverlib_cell2d_impl.h \ + private/iricsolverlib_grid2d_impl.h \ + private/iricsolverlib_point2d_detail.h SOURCES += filelocker.cpp \ iric_ftoc.c \ iriclib.cpp \ @@ -75,6 +87,13 @@ SOURCES += filelocker.cpp \ iriclib_polyline.cpp \ iriclib_riversurvey.cpp \ iriclib_single.c \ + iricsolver_ftoc.c \ + iricsolverlib.cpp \ + iricsolverlib_cell2d.cpp \ + iricsolverlib_grid2d.cpp \ + iricsolverlib_quadcell.cpp \ + iricsolverlib_rect2d.cpp \ + iricsolverlib_tricell.cpp \ private/iriclib_cgnsfile_solutionwriter.cpp \ private/iriclib_cgnsfile_solutionwriterdividesolutions.cpp \ private/iriclib_cgnsfile_solutionwriterstandard.cpp diff --git a/iricsolver_ftoc.c b/iricsolver_ftoc.c new file mode 100644 index 0000000..af3e846 --- /dev/null +++ b/iricsolver_ftoc.c @@ -0,0 +1,105 @@ +#include +#include + +#include "iricsolverlib.h" +#include "iricsolverlib_api.h" + +#include "fortran_macros.h" + +#include + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *\ + * Convert between Fortran and C strings * +\* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +static void string_2_C_string(char *string, int string_length, + char *c_string, int max_len, int *ierr) { + int i, iend; + + if (string == NULL || c_string == NULL) { +// cgi_error ("NULL string pointer"); + *ierr = CG_ERROR; + return; + } + + /** Skip and trailing blanks **/ + for (iend = string_length-1; iend >= 0; iend--) { + if (string[iend] != ' ') break; + } + if (iend >= max_len) iend = max_len - 1; + + /** Copy the non-trailing blank portion of the string **/ + for (i = 0; i <= iend; i++) + c_string[i] = string[i]; + + /** NULL terminate the C string **/ + c_string[i] = '\0'; + *ierr = CG_OK; +} + +static void string_2_F_string(char *c_string, char *string, + int string_length, int *ierr) { + int i; + size_t len; + + if (c_string == NULL || string == NULL) { +// cgi_error ("NULL string pointer"); + *ierr = CG_ERROR; + return; + } + len = strlen(c_string); + if (len > string_length) len = string_length; + + for (i = 0; i < len; i++) + string[i] = c_string[i]; + while (i < string_length) + string[i++] = ' '; + *ierr = CG_OK; +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_open_f, IRIC_SOLVER_GRID2D_OPEN_F) (int *fin, int *baseId, int *zoneId, int *gridId, int* handle, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Open(*fin, *baseId, *zoneId, *gridId, handle); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_read_cellcount_f, IRIC_SOLVER_GRID2D_READ_CELLCOUNT_F) (int *gridHandle, int *cellCount, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Read_CellCount(*gridHandle, cellCount); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_read_cellnodecount_f, IRIC_SOLVER_GRID2D_READ_CELLNODECOUNT_F) (int *gridHandle, int *cellId, int* nodeCount, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Read_CellNodeCount(*gridHandle, *cellId, nodeCount); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_read_cellnodeids_f, IRIC_SOLVER_GRID2D_READ_CELLNODEIDS_F) (int *gridHandle, int *cellId, int* nodeIds, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Read_CellNodeIds(*gridHandle, *cellId, nodeIds); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_read_cellarea_f, IRIC_SOLVER_GRID2D_READ_CELLAREA_F) (int *gridHandle, int *cellId, double* area, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Read_CellArea(*gridHandle, *cellId, area); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_getregion_f, IRIC_SOLVER_GRID2D_GETREGION_F) (int *gridHandle, double *xmin, double *xmax, double *ymin, double *ymax, int *ier) +{ + *ier = iRIC_Solver_Grid2D_GetRegion(*gridHandle, xmin, xmax, ymin, ymax); +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_interpolate_f, IRIC_SOLVER_GRID2D_INTERPOLATE_F) (int *gridHandle, double *x, double *y, int *ok, int *count, int* nodeids, double* weights, int *ier) +{ + int i; + size_t tmpNodeIds[4]; + + *ier = iRIC_Solver_Grid2D_Interpolate(*gridHandle, *x, *y, ok, count, &(tmpNodeIds[0]), weights); + + for (i = 0; i < *count; ++i) { + *(nodeids + i) = (int) (tmpNodeIds[i]); + } +} + +void IRICSOLVERLIB_API FMNAME(iric_solver_grid2d_close_f, IRIC_SOLVER_GRID2D_CLOSE_F) (int *gridHandle, int *ier) +{ + *ier = iRIC_Solver_Grid2D_Close(*gridHandle); +} diff --git a/iricsolverlib.cpp b/iricsolverlib.cpp new file mode 100644 index 0000000..79b1f70 --- /dev/null +++ b/iricsolverlib.cpp @@ -0,0 +1,125 @@ +#include +#include "iricsolverlib_cell2d.h" +#include "iricsolverlib_grid2d.h" +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_rect2d.h" + +static std::vector grid2ds; + +using namespace iRICSolverLib; + +extern "C" { + +int iRIC_Solver_Grid2D_Open(int fin, int baseId, int zoneId, int gridId, int* handle) +{ + Grid2D* grid = new Grid2D(); + grid->load(fin, baseId, zoneId, gridId); + grid2ds.push_back(grid); + *handle = static_cast (grid2ds.size()); + + return 0; +} + +int iRIC_Solver_Grid2D_Read_CellCount(int gridHandle, int* cellCount) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + *cellCount = grid->cellCount(); + + return 0; +} + +int iRIC_Solver_Grid2D_Read_CellNodeCount(int gridHandle, int cellId, int* nodeCount) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + if (cellId < 1) {return -3;} + if (cellId > grid->cellCount()) {return -3;} + + Cell2D* cell = grid->cell(static_cast (cellId)); + *nodeCount = cell->nodeCount(); + + return 0; +} + +int iRIC_Solver_Grid2D_Read_CellNodeIds(int gridHandle, int cellId, int* nodeIds) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + if (cellId < 1) {return -3;} + if (cellId > grid->cellCount()) {return -3;} + + Cell2D* cell = grid->cell(static_cast (cellId)); + for (int i = 0; i < cell->nodeCount(); ++i) { + *(nodeIds + i) = static_cast (cell->nodeId(i + 1)); + } + + return 0; +} + +int iRIC_Solver_Grid2D_Read_CellArea(int gridHandle, int cellId, double* area) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + if (cellId < 1) {return -3;} + if (cellId > grid->cellCount()) {return -3;} + + Cell2D* cell = grid->cell(static_cast (cellId)); + *area = cell->area(); + + return 0; +} + +int iRIC_Solver_Grid2D_GetRegion(int gridHandle, double* xmin, double* xmax, double* ymin, double* ymax) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + Rect2D rect = grid->boundingRect(); + *xmin = rect.xMin(); + *xmax = rect.xMax(); + *ymin = rect.yMin(); + *ymax = rect.yMax(); + + return 0; +} + +int iRIC_Solver_Grid2D_Interpolate(int gridHandle, double x, double y, int* ok, int* count, size_t* nodeids, double* weights) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + Grid2D* grid = grid2ds.at(gridHandle - 1); + if (grid == 0) {return -2;} + + *ok = 0; + bool success = grid->interpolate(Point2D(x, y), count, nodeids, weights); + if (success) { + *ok = 1; + } + return 0; +} + +int iRIC_Solver_Grid2D_Close(int gridHandle) +{ + if (static_cast (gridHandle) > grid2ds.size()) {return -1;} + + delete grid2ds[gridHandle - 1]; + grid2ds[gridHandle - 1] = 0; + return 0; +} + +} // extern "C" diff --git a/iricsolverlib.h b/iricsolverlib.h new file mode 100644 index 0000000..81361b5 --- /dev/null +++ b/iricsolverlib.h @@ -0,0 +1,22 @@ +#include "iricsolverlib_api.h" + +#ifdef __cplusplus +extern "C" { +#endif + +// 2D Grid related functions +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Open(int fin, int baseId, int zoneId, int gridId, int* handle); + +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Read_CellCount(int gridHandle, int* cellCount); +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Read_CellNodeCount(int gridHandle, int cellId, int* nodeCount); +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Read_CellNodeIds(int gridHandle, int cellId, int* nodeIds); +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Read_CellArea(int gridHandle, int cellId, double* area); + +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_GetRegion(int gridHandle, double* xmin, double* xmax, double* ymin, double* ymax); +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Interpolate(int gridHandle, double x, double y, int* ok, int* count, size_t* nodeids, double* weights); + +int IRICSOLVERLIB_API iRIC_Solver_Grid2D_Close(int gridHandle); + +#ifdef __cplusplus +} +#endif diff --git a/iricsolverlib_api.h b/iricsolverlib_api.h new file mode 100644 index 0000000..94c2e63 --- /dev/null +++ b/iricsolverlib_api.h @@ -0,0 +1,16 @@ +#ifndef IRICSOLVERLIB_API +# ifdef _WIN32 +# if defined(IRICSOLVERLIBSTATIC_LIBRARY) +# define IRICSOLVERLIB_API +# else +# if defined(IRICSOLVERLIBDLL_LIBRARY) +# define IRICSOLVERLIB_API __declspec(dllexport) +# else +# define IRICSOLVERLIB_API __declspec(dllimport) +# endif +# endif +# else +# define IRICSOLVERLIB_API +# endif +#endif + diff --git a/iricsolverlib_cell2d.cpp b/iricsolverlib_cell2d.cpp new file mode 100644 index 0000000..8e6b67a --- /dev/null +++ b/iricsolverlib_cell2d.cpp @@ -0,0 +1,55 @@ +#include "iricsolverlib_cell2d.h" +#include "iricsolverlib_grid2d.h" +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_rect2d.h" +#include "private/iricsolverlib_cell2d_impl.h" + +using namespace iRICSolverLib; + +Cell2D::Impl::Impl(Grid2D* const grid) : + m_grid (grid) +{} + +// public interfaces + +Cell2D::Cell2D(Grid2D* const grid) : + impl (new Impl(grid)) +{} + +Cell2D::~Cell2D() +{ + delete impl; +} + +int Cell2D::nodeCount() const +{ + return static_cast (impl->m_nodeIds.size()); +} + +size_t Cell2D::nodeId(int id) const +{ + return impl->m_nodeIds.at(id - 1) + 1; +} + +Point2D Cell2D::node(int id) const +{ + return impl->m_grid->node(nodeId(id)); +} + +Rect2D Cell2D::boundingRect() const +{ + double xmin, xmax, ymin, ymax; + for (size_t i = 0; i < nodeCount(); ++i) { + Point2D p = node(i + 1); + if (i == 0 || xmin > p.x()) {xmin = p.x();} + if (i == 0 || xmax < p.x()) {xmax = p.x();} + if (i == 0 || ymin > p.y()) {ymin = p.y();} + if (i == 0 || ymax < p.y()) {ymax = p.y();} + } + return Rect2D(xmin, xmax, ymin, ymax); +} + +void Cell2D::addNode(size_t id) +{ + impl->m_nodeIds.push_back(id - 1); +} diff --git a/iricsolverlib_cell2d.h b/iricsolverlib_cell2d.h new file mode 100644 index 0000000..7d7c5ce --- /dev/null +++ b/iricsolverlib_cell2d.h @@ -0,0 +1,47 @@ +#ifndef IRICSOLVERLIB_CELL2D_H +#define IRICSOLVERLIB_CELL2D_H + +#include "iricsolverlib_api.h" + +#include + +namespace iRICSolverLib { + +class Grid2D; +class Point2D; +class Rect2D; + +class IRICSOLVERLIB_API Cell2D +{ +public: + Cell2D(Grid2D* const grid); + virtual ~Cell2D(); + + int nodeCount() const; + size_t nodeId(int id) const; + Point2D node(int id) const; + + Rect2D boundingRect() const; + + virtual bool interpolate(const Point2D& point, double* weight) const = 0; + virtual double area() const = 0; + +protected: + void addNode(size_t id); + +private: + class Impl; + Impl* impl; + + // copy not allowed + Cell2D(const Cell2D& cell); + Cell2D& operator=(const Cell2D& cell); +}; + +} // iRICSolverLib + +#ifdef _DEBUG + #include "private/iricsolverlib_cell2d_impl.h" +#endif // _DEBUG + +#endif // IRICSOLVERLIB_CELL2D_H diff --git a/iricsolverlib_grid2d.cpp b/iricsolverlib_grid2d.cpp new file mode 100644 index 0000000..952fc6e --- /dev/null +++ b/iricsolverlib_grid2d.cpp @@ -0,0 +1,311 @@ +#include "iricsolverlib_cell2d.h" +#include "iricsolverlib_grid2d.h" +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_quadcell.h" +#include "iricsolverlib_rect2d.h" +#include "iricsolverlib_tricell.h" +#include "private/iricsolverlib_grid2d_impl.h" + +#include +#include +#include + +#include + +using namespace iRICSolverLib; + +namespace { + +} // namespace + +Grid2D::Impl::Impl(Grid2D* grid) : + m_grid (grid) +{} + +Grid2D::Impl::~Impl() +{ + for (size_t i = 0; i < m_cells.size(); ++i) { + delete m_cells.at(i); + } +} + +void Grid2D::Impl::loadStructuredGrid(int fn, int baseId, int zoneId, int gridId) +{ + char zoneName[32]; + int size[9]; + int ier; + + cg_zone_read(fn, baseId, zoneId, zoneName, &(size[0])); + + int nodeCount = size[0] * size[1]; + std::vector xVec(nodeCount, 0); + std::vector yVec(nodeCount, 0); + + ier = cg_goto(fn, baseId, "Zone_t", zoneId, "GridCoordinates_t", gridId, "end"); + // load X and Y + cg_array_read_as(1, RealDouble, &(xVec.front())); + cg_array_read_as(2, RealDouble, &(yVec.front())); + + m_nodes.clear(); + m_nodes.reserve(nodeCount); + for (int i = 0; i < nodeCount; ++i) { + m_nodes.push_back(Point2D(xVec.at(i), yVec.at(i))); + } + + for (int i = 0; i < size[0] - 1; ++i) { + for (int j = 0; j < size[1] - 1; ++j) { + size_t id0 = structuredIndex(i , j , size); + size_t id1 = structuredIndex(i + 1, j , size); + size_t id2 = structuredIndex(i + 1, j + 1, size); + size_t id3 = structuredIndex(i , j + 1, size); + m_cells.push_back(new QuadCell(id0 + 1, id1 + 1, id2 + 1, id3 + 1, m_grid)); + } + } +} + +void Grid2D::Impl::insertTriangleCells(int fn, int baseId, int zoneId, int secId) +{ + int numCells; + cg_ElementDataSize(fn, baseId, zoneId, secId, &numCells); + numCells = numCells / 3; + std::vector elements(3 * numCells, 0); + cg_elements_read(fn, baseId, zoneId, secId, &(elements.front()), NULL); + + for (int i = 0; i < numCells; ++i) { + size_t id0 = elements.at(i * 3); + size_t id1 = elements.at(i * 3 + 1); + size_t id2 = elements.at(i * 3 + 2); + m_cells.push_back(new TriCell(id0, id1, id2, m_grid)); + } +} + +void Grid2D::Impl::insertQuadCells(int fn, int baseId, int zoneId, int secId) +{ + int numCells; + cg_ElementDataSize(fn, baseId, zoneId, secId, &numCells); + numCells = numCells / 4; + std::vector elements(4 * numCells, 0); + cg_elements_read(fn, baseId, zoneId, secId, &(elements.front()), NULL); + + for (int i = 0; i < numCells; ++i) { + size_t id0 = elements.at(i * 4); + size_t id1 = elements.at(i * 4 + 1); + size_t id2 = elements.at(i * 4 + 2); + size_t id3 = elements.at(i * 4 + 3); + m_cells.push_back(new QuadCell(id0, id1, id2, id3, m_grid)); + } +} + +void Grid2D::Impl::loadUnstructuredGrid(int fn, int baseId, int zoneId, int gridId) +{ + char zoneName[32]; + int size[9]; + int ier; + + cg_zone_read(fn, baseId, zoneId, zoneName, &(size[0])); + + int nodeCount = size[0]; + std::vector xVec(nodeCount, 0); + std::vector yVec(nodeCount, 0); + + ier = cg_goto(fn, baseId, "Zone_t", zoneId, "GridCoordinates_t", gridId, "end"); + // load X and Y + cg_array_read_as(1, RealDouble, &(xVec.front())); + cg_array_read_as(2, RealDouble, &(yVec.front())); + + m_nodes.clear(); + m_nodes.reserve(nodeCount); + for (int i = 0; i < nodeCount; ++i) { + m_nodes.push_back(Point2D(xVec.at(i), yVec.at(i))); + } + + // load sections + int numSections; + cg_nsections(fn, baseId, zoneId, &numSections); + for (int S = 1; S <= numSections; ++S) { + ElementType_t eType; + int startIndex, endIndex; + int nBndry, parent_flag; + char buffer[32]; + cg_section_read(fn, baseId, zoneId, S, buffer, &eType, &startIndex, &endIndex, &nBndry, &parent_flag); + if (eType == TRI_3) { + // Triangle + insertTriangleCells(fn, baseId, zoneId, S); + } else if (eType == QUAD_4) { + // Quadrangle + insertQuadCells(fn, baseId, zoneId, S); + } + } +} + +void Grid2D::Impl::setupBackGrid() +{ + m_backGridX.clear(); + m_backGridY.clear(); + m_backGridCells.clear(); + + Rect2D rect = m_grid->boundingRect(); + size_t numCells = m_grid->cellCount(); + int divNum = static_cast (std::sqrt(static_cast(numCells) / 4)); + if (divNum == 0) {divNum = 1;} + + double xDelta = (rect.xMax() - rect.xMin()) / divNum; + double yDelta = (rect.yMax() - rect.yMin()) / divNum; + for (int i = 0; i < divNum; ++i) { + m_backGridX.push_back(rect.xMin() + i * xDelta); + m_backGridY.push_back(rect.yMin() + i * yDelta); + } + m_backGridX.push_back(rect.xMax()); + m_backGridY.push_back(rect.yMax()); + + std::vector emptyVec; + for (int i = 0; i < divNum; ++i) { + for (int j = 0; j < divNum; ++j) { + m_backGridCells.push_back(emptyVec); + } + } + + for (size_t i = 0; i < m_cells.size(); ++i) { + Cell2D* cell = m_cells.at(i); + Rect2D rect = cell->boundingRect(); + + std::vector::const_iterator x_lb = std::lower_bound(m_backGridX.begin(), m_backGridX.end(), rect.xMin()); + if (*x_lb != rect.xMin() && x_lb != m_backGridX.begin()) {--x_lb;} + + std::vector::const_iterator x_ub = std::lower_bound(m_backGridX.begin(), m_backGridX.end(), rect.xMax()); + + std::vector::const_iterator y_lb = std::lower_bound(m_backGridY.begin(), m_backGridY.end(), rect.yMin()); + if (*y_lb != rect.yMin() && y_lb != m_backGridY.begin()) {--y_lb;} + + std::vector::const_iterator y_ub = std::lower_bound(m_backGridY.begin(), m_backGridY.end(), rect.yMax()); + + for (std::vector::const_iterator xit = x_lb; xit != x_ub; ++xit) { + size_t x_idx = xit - m_backGridX.begin(); + for (std::vector::const_iterator yit = y_lb; yit != y_ub; ++yit) { + size_t y_idx = yit - m_backGridY.begin(); + + size_t idx = x_idx + y_idx * (m_backGridX.size() - 1); + m_backGridCells[idx].push_back(cell); + } + } + } +} + +size_t Grid2D::Impl::structuredIndex(size_t i, size_t j, int* size) const +{ + return i + j * (*size); +} + +// ----------------------------------------------------------- +// public interfaces +// ----------------------------------------------------------- + +Grid2D::Grid2D() +{ + impl = new Impl(this); +} + +Grid2D::~Grid2D() +{ + delete impl; +} + +void Grid2D::load(int cgnsIn, int baseId, int zoneId, int gridId) +{ + ZoneType_t zType; + + cg_zone_type(cgnsIn, baseId, zoneId, &zType); + if (zType == Structured) { + impl->loadStructuredGrid(cgnsIn, baseId, zoneId, gridId); + } else { + impl->loadUnstructuredGrid(cgnsIn, baseId, zoneId, gridId); + } + impl->setupBackGrid(); +} + +size_t Grid2D::nodeCount() const +{ + return impl->m_nodes.size(); +} + +Point2D Grid2D::node(size_t nodeId) const +{ + return impl->m_nodes.at(nodeId - 1); +} + +size_t Grid2D::cellCount() const +{ + return impl->m_cells.size(); +} + +Cell2D* Grid2D::cell(size_t cellId) const +{ + return impl->m_cells.at(cellId - 1); +} + +Rect2D Grid2D::boundingRect() const +{ + Rect2D ret; + for (int i = 0; i < impl->m_cells.size(); ++i) { + Rect2D cellRect = impl->m_cells.at(i)->boundingRect(); + if (i == 0) { + ret = cellRect; + } else { + ret = ret + cellRect; + } + } + return ret; +} + +bool Grid2D::interpolate(const Point2D& point, int *count, size_t* nodeIds, double* weight) const +{ + std::vector::iterator lb_x = std::lower_bound(impl->m_backGridX.begin(), impl->m_backGridX.end(), point.x()); + if (lb_x == impl->m_backGridX.end()) {return false;} + if (*lb_x != point.x() && lb_x != impl->m_backGridX.begin()) {-- lb_x;} + + std::vector::iterator lb_y = std::lower_bound(impl->m_backGridY.begin(), impl->m_backGridY.end(), point.y()); + if (lb_y == impl->m_backGridY.end()) {return false;} + if (*lb_y != point.y() && lb_y != impl->m_backGridY.begin()) {--lb_y;} + + size_t bgCellIdx = (lb_x - impl->m_backGridX.begin()) + + (lb_y - impl->m_backGridY.begin()) * (impl->m_backGridX.size() - 1); + + const std::vector& cells = impl->m_backGridCells.at(bgCellIdx); + for (int i = 0; i < cells.size(); ++i) { + Cell2D* cell = cells.at(i); + *count = cell->nodeCount(); + bool ok = cell->interpolate(point, weight); + if (ok) { + for (int j = 0; j < cell->nodeCount(); ++j) { + *(nodeIds + j) = cell->nodeId(j + 1); + } + return true; + } + } + return false; +} + +void Grid2D::addNode(const Point2D& p) +{ + impl->m_nodes.push_back(p); +} + +void Grid2D::addNode(double x, double y) +{ + addNode(Point2D(x, y)); +} + +void Grid2D::addTriCell(size_t id1, size_t id2, size_t id3) +{ + impl->m_cells.push_back(new TriCell(id1, id2, id3, this)); +} + +void Grid2D::addQuadCell(size_t id1, size_t id2, size_t id3, size_t id4) +{ + impl->m_cells.push_back(new QuadCell(id1, id2, id3, id4, this)); +} + +void Grid2D::setupBackGrid() +{ + impl->setupBackGrid(); +} diff --git a/iricsolverlib_grid2d.h b/iricsolverlib_grid2d.h new file mode 100644 index 0000000..4fb4ea0 --- /dev/null +++ b/iricsolverlib_grid2d.h @@ -0,0 +1,51 @@ +#ifndef IRICSOLVERLIB_GRID2D_H +#define IRICSOLVERLIB_GRID2D_H + +#include "iricsolverlib_api.h" + +#include + +namespace iRICSolverLib { + +class Cell2D; +class Point2D; +class Rect2D; + +class IRICSOLVERLIB_API Grid2D +{ +public: + Grid2D(); + ~Grid2D(); + + void load(int cgnsIn, int baseId, int zoneId, int gridId); + + size_t nodeCount() const; + Point2D node(size_t nodeId) const; + + size_t cellCount() const; + Cell2D* cell(size_t cellId) const; + + Rect2D boundingRect() const; + + bool interpolate(const Point2D& point, int *count, size_t* nodeIds, double* weight) const; + + void addNode(const Point2D& p); + void addNode(double x, double y); + + void addTriCell(size_t id1, size_t id2, size_t id3); + void addQuadCell(size_t id1, size_t id2, size_t id3, size_t id4); + + void setupBackGrid(); + +private: + class Impl; + Impl* impl; +}; + +} // iRICSolverLib + +#ifdef _DEBUG + #include "private/iricsolverlib_grid2d_impl.h" +#endif // _DEBUG + +#endif // IRICSOLVERLIB_GRID2D_H diff --git a/iricsolverlib_point2d.h b/iricsolverlib_point2d.h new file mode 100644 index 0000000..60caaf7 --- /dev/null +++ b/iricsolverlib_point2d.h @@ -0,0 +1,30 @@ +#ifndef IRICSOLVERLIB_POINT2D_H +#define IRICSOLVERLIB_POINT2D_H + +#include "iricsolverlib_api.h" + +namespace iRICSolverLib { + +class Point2D +{ +public: + Point2D(); + Point2D(double x, double y); + + double x() const; + double y() const; + +private: + double m_x; + double m_y; +}; + +Point2D operator+(const Point2D& p1, const Point2D& p2); +Point2D operator-(const Point2D& p1, const Point2D& p2); +double operator*(const Point2D& p1, const Point2D& p2); + +} // iRICSolverLib + +#include "private/iricsolverlib_point2d_detail.h" + +#endif // IRICSOLVERLIB_POINT2D_H diff --git a/iricsolverlib_quadcell.cpp b/iricsolverlib_quadcell.cpp new file mode 100644 index 0000000..f6837b1 --- /dev/null +++ b/iricsolverlib_quadcell.cpp @@ -0,0 +1,75 @@ +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_quadcell.h" +#include "iricsolverlib_rect2d.h" +#include "iricsolverlib_tricell.h" + +#include + +using namespace iRICSolverLib; + +namespace { + +const double INSIDE_DELTA = 1.0E-8; + +bool interpolateTri(const QuadCell& cell, const Point2D& point, int id1, int id2, int id3, double* weights) +{ + double s, t, u; + TriCell::calcSTU(point, cell.node(id1), cell.node(id2), cell.node(id3), &s, &t, &u); + + if (s < - INSIDE_DELTA || t < - INSIDE_DELTA || u < - INSIDE_DELTA) { + return false; + } + + *(weights) = s; + *(weights + 1) = t; + *(weights + 2) = u; + return true; +} + +} // namespace + +QuadCell::QuadCell(size_t id1, size_t id2, size_t id3, size_t id4, Grid2D* const grid) : + Cell2D(grid) +{ + addNode(id1); + addNode(id2); + addNode(id3); + addNode(id4); +} + +QuadCell::~QuadCell() +{} + +bool QuadCell::interpolate(const Point2D& point, double* weight) const +{ + if (! boundingRect().contains(point)) {return false;} + + double tmpWeights[3]; + bool ok; + // try triangle with, 0, 1, 2 + ok = interpolateTri(*this, point, 1, 2, 3, &(tmpWeights[0])); + if (ok) { + * (weight ) = tmpWeights[0]; + * (weight + 1) = tmpWeights[1]; + * (weight + 2) = tmpWeights[2]; + * (weight + 3) = 0; + return true; + } + + // try triangle with 0, 2, 3 + ok = interpolateTri(*this, point, 1, 3, 4, &(tmpWeights[0])); + if (ok) { + * (weight ) = tmpWeights[0]; + * (weight + 1) = 0; + * (weight + 2) = tmpWeights[1]; + * (weight + 3) = tmpWeights[2]; + return true; + } + return false; +} + +double QuadCell::area() const +{ + return TriCell::calcArea(node(1), node(2), node(3)) + + TriCell::calcArea(node(1), node(3), node(4)); +} diff --git a/iricsolverlib_quadcell.h b/iricsolverlib_quadcell.h new file mode 100644 index 0000000..0eb1af6 --- /dev/null +++ b/iricsolverlib_quadcell.h @@ -0,0 +1,21 @@ +#ifndef IRICSOLVERLIB_QUADCELL_H +#define IRICSOLVERLIB_QUADCELL_H + +#include "iricsolverlib_api.h" +#include "iricsolverlib_cell2d.h" + +namespace iRICSolverLib { + +class IRICSOLVERLIB_API QuadCell : public Cell2D +{ +public: + QuadCell(size_t id1, size_t id2, size_t id3, size_t id4, Grid2D* const grid); + ~QuadCell(); + + bool interpolate(const Point2D& point, double* weight) const; + double area() const; +}; + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_QUADCELL_H diff --git a/iricsolverlib_rect2d.cpp b/iricsolverlib_rect2d.cpp new file mode 100644 index 0000000..62d6a41 --- /dev/null +++ b/iricsolverlib_rect2d.cpp @@ -0,0 +1,120 @@ +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_rect2d.h" + +#include + +using namespace std; + +namespace iRICSolverLib { + +Rect2D::Rect2D() : + m_xMin (0), + m_xMax (0), + m_yMin (0), + m_yMax (0) +{} + +Rect2D::Rect2D(double xmin, double xmax, double ymin, double ymax) : + m_xMin (xmin), + m_xMax (xmax), + m_yMin (ymin), + m_yMax (ymax) +{} + +Rect2D::Rect2D(const Rect2D& rect) : + m_xMin (rect.xMin()), + m_xMax (rect.xMax()), + m_yMin (rect.yMin()), + m_yMax (rect.yMax()) +{} + +Rect2D& Rect2D::operator=(const Rect2D& rect) +{ + m_xMin = rect.xMin(); + m_xMax = rect.xMax(); + m_yMin = rect.yMin(); + m_yMax = rect.yMax(); + + return *this; +} + +double Rect2D::xMin() const +{ + return m_xMin; +} + +double Rect2D::xMax() const +{ + return m_xMax; +} + +double Rect2D::yMin() const +{ + return m_yMin; +} + +double Rect2D::yMax() const +{ + return m_yMax; +} + +double Rect2D::xWidth() const +{ + return m_xMax - m_xMin; +} + +double Rect2D::yWidth() const +{ + return m_yMax - m_yMin; +} + +bool Rect2D::contains(double x, double y) const +{ + return (m_xMin <= x && m_xMax >= x && m_yMin <= y && m_yMax >= y); +} + +bool Rect2D::contains(const Point2D& point) const +{ + return contains(point.x(), point.y()); +} + +bool Rect2D::intersects(const Rect2D& rect) const +{ + if (m_xMin > rect.xMax()) {return false;} + if (m_xMax < rect.xMin()) {return false;} + if (m_yMin > rect.yMax()) {return false;} + if (m_yMax < rect.yMin()) {return false;} + + return true; +} + +void Rect2D::setXMin(double xmin) +{ + m_xMin = xmin; +} + +void Rect2D::setXMax(double xmax) +{ + m_xMax = xmax; +} + +void Rect2D::setYMin(double ymin) +{ + m_yMin = ymin; +} + +void Rect2D::setYMax(double ymax) +{ + m_yMax = ymax; +} + +Rect2D operator+(const Rect2D& rect1, const Rect2D& rect2) +{ + return Rect2D( + min(rect1.xMin(), rect2.xMin()), + max(rect1.xMax(), rect2.xMax()), + min(rect1.yMin(), rect2.yMin()), + max(rect1.yMax(), rect2.yMax())); +} + +} // namespace iRICSolverLib diff --git a/iricsolverlib_rect2d.h b/iricsolverlib_rect2d.h new file mode 100644 index 0000000..a867549 --- /dev/null +++ b/iricsolverlib_rect2d.h @@ -0,0 +1,47 @@ +#ifndef IRICSOLVERLIB_RECT2D_H +#define IRICSOLVERLIB_RECT2D_H + +#include "iricsolverlib_api.h" + +namespace iRICSolverLib { + +class Point2D; + +class IRICSOLVERLIB_API Rect2D +{ +public: + Rect2D(); + Rect2D(double xmin, double xmax, double ymin, double ymax); + Rect2D(const Rect2D& rect); + + Rect2D& operator=(const Rect2D& rect); + + double xMin() const; + double xMax() const; + double yMin() const; + double yMax() const; + + double xWidth() const; + double yWidth() const; + + bool contains(double x, double y) const; + bool contains(const Point2D& point) const; + bool intersects(const Rect2D& rect) const; + + void setXMin(double xmin); + void setXMax(double xmax); + void setYMin(double ymin); + void setYMax(double ymax); + +private: + double m_xMin; + double m_xMax; + double m_yMin; + double m_yMax; +}; + +Rect2D operator +(const Rect2D& rect1, const Rect2D& rect2); + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_RECT2D_H diff --git a/iricsolverlib_tricell.cpp b/iricsolverlib_tricell.cpp new file mode 100644 index 0000000..b391bf2 --- /dev/null +++ b/iricsolverlib_tricell.cpp @@ -0,0 +1,64 @@ +#include "iricsolverlib_tricell.h" +#include "iricsolverlib_point2d.h" +#include "iricsolverlib_rect2d.h" + +#include +#include + +using namespace iRICSolverLib; + +namespace { + +const double INSIDE_DELTA = 1.0E-8; + +} // namespace + +TriCell::TriCell(size_t id1, size_t id2, size_t id3, Grid2D* const grid) : + Cell2D(grid) +{ + addNode(id1); + addNode(id2); + addNode(id3); +} + +TriCell::~TriCell() +{} + +bool TriCell::interpolate(const Point2D& point, double* weight) const +{ + if (! boundingRect().contains(point)) {return false;} + + double s, t, u; + calcSTU(point, node(1), node(2), node(3), &s, &t, &u); + + if (s < - INSIDE_DELTA || t < - INSIDE_DELTA || u < - INSIDE_DELTA) { + return false; + } + + *(weight) = s; + *(weight + 1) = t; + *(weight + 2) = u; + return true; +} + +double TriCell::area() const +{ + return calcArea(node(1), node(2), node(3)); +} + +void TriCell::calcSTU(const Point2D& point, const Point2D& node0, const Point2D& node1,const Point2D& node2, double* s, double* t, double* u) +{ + Point2D p = point - node0; + Point2D p1 = node1 - node0; + Point2D p2 = node2 - node0; + + double m = p1.x() * p2.y() - p1.y() * p2.x(); + *t = (p.x() * p2.y() - p.y() * p2.x()) / m; + *u = - (p.x() * p1.y() - p.y() * p1.x()) / m; + *s = (1 - *t - *u); +} + +double TriCell::calcArea(const Point2D& p1, const Point2D& p2, const Point2D& p3) +{ + return 0.5 * std::fabs((p2.x() - p1.x()) * (p3.y() - p1.y()) - (p2.y() - p1.y()) * (p3.x() - p1.x())); +} diff --git a/iricsolverlib_tricell.h b/iricsolverlib_tricell.h new file mode 100644 index 0000000..102400b --- /dev/null +++ b/iricsolverlib_tricell.h @@ -0,0 +1,24 @@ +#ifndef IRICSOLVERLIB_TRICELL_H +#define IRICSOLVERLIB_TRICELL_H + +#include "iricsolverlib_api.h" +#include "iricsolverlib_cell2d.h" + +namespace iRICSolverLib { + +class IRICSOLVERLIB_API TriCell : public Cell2D +{ +public: + TriCell(size_t id1, size_t id2, size_t id3, Grid2D* const grid); + ~TriCell(); + + bool interpolate(const Point2D& point, double* weight) const; + double area() const; + + static void calcSTU(const Point2D& point, const Point2D& p1, const Point2D& p2,const Point2D& p3, double* s, double* t, double* u); + static double calcArea(const Point2D& p1, const Point2D& p2, const Point2D& p3); +}; + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_TRICELL_H diff --git a/private/iricsolverlib_cell2d_impl.h b/private/iricsolverlib_cell2d_impl.h new file mode 100644 index 0000000..adbf445 --- /dev/null +++ b/private/iricsolverlib_cell2d_impl.h @@ -0,0 +1,22 @@ +#ifndef IRICSOLVERLIB_CELL2D_IMPL_H +#define IRICSOLVERLIB_CELL2D_IMPL_H + +#include "../iricsolverlib_cell2d.h" + +#include + +namespace iRICSolverLib { + +class Cell2D::Impl +{ +public: + Impl(Grid2D* const grid); + + Grid2D* const m_grid; + std::vector m_nodeIds; +}; + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_CELL2D_IMPL_H + diff --git a/private/iricsolverlib_grid2d_impl.h b/private/iricsolverlib_grid2d_impl.h new file mode 100644 index 0000000..071574a --- /dev/null +++ b/private/iricsolverlib_grid2d_impl.h @@ -0,0 +1,38 @@ +#ifndef IRICSOLVERLIB_GRID2D_IMPL_H +#define IRICSOLVERLIB_GRID2D_IMPL_H + +#include "../iricsolverlib_grid2d.h" + +#include + +namespace iRICSolverLib { + +class Grid2D::Impl +{ +public: + Impl(Grid2D* grid); + ~Impl(); + void loadStructuredGrid(int fn, int baseId, int zoneId, int gridId); + void loadUnstructuredGrid(int fn, int baseId, int zoneId, int gridId); + void setupBackGrid(); + + size_t structuredIndex(size_t i, size_t j, int* size) const; + +private: + void insertTriangleCells(int fn, int baseId, int zoneId, int secId); + void insertQuadCells(int fn, int baseId, int zoneId, int secId); + +public: + Grid2D* m_grid; + std::vector m_nodes; + std::vector m_cells; + + std::vector m_backGridX; + std::vector m_backGridY; + + std::vector< std::vector > m_backGridCells; +}; + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_GRID2D_IMPL_H diff --git a/private/iricsolverlib_point2d_detail.h b/private/iricsolverlib_point2d_detail.h new file mode 100644 index 0000000..e3c9891 --- /dev/null +++ b/private/iricsolverlib_point2d_detail.h @@ -0,0 +1,43 @@ +#ifndef IRICSOLVERLIB_POINT2D_DETAIL_H +#define IRICSOLVERLIB_POINT2D_DETAIL_H + +#include "../iricsolverlib_point2d.h" + +namespace iRICSolverLib { + +inline Point2D::Point2D() : + m_x (0), m_y(0) +{} + +inline Point2D::Point2D(double x, double y) : + m_x (x), m_y(y) +{} + +inline double Point2D::x() const +{ + return m_x; +} + +inline double Point2D::y() const +{ + return m_y; +} + +inline Point2D operator+(const Point2D& p1, const Point2D& p2) +{ + return Point2D(p1.x() + p2.x(), p1.y() + p2.y()); +} + +inline Point2D operator-(const Point2D& p1, const Point2D& p2) +{ + return Point2D(p1.x() - p2.x(), p1.y() - p2.y()); +} + +inline double operator*(const Point2D& p1, const Point2D& p2) +{ + return p1.x() * p2.y() - p1.y() * p2.x(); +} + +} // iRICSolverLib + +#endif // IRICSOLVERLIB_POINT2D_DETAIL_H