diff --git a/CMakeLists.txt b/CMakeLists.txt index c0ebb2b..f4e539e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,7 @@ endif () set(WONTON_VERSION_MAJOR 1) set(WONTON_VERSION_MINOR 2) -set(WONTON_VERSION_PATCH 1) +set(WONTON_VERSION_PATCH 2) # Top level target @@ -95,8 +95,8 @@ if (WONTON_ENABLE_FleCSI AND NOT FleCSI_LIBRARIES) message(STATUS "FleCSI_LIBRARIES=${FleCSI_LIBRARIES}" ) message(STATUS "FleCSISP_LIBRARIES=${FleCSISP_LIBRARIES}" ) - target_include_directories(wonton INTERFACE ${FleCSI_INCLUDE_DIR}) - target_include_directories(wonton INTERFACE ${FleCSISP_INCLUDE_DIR}) + target_include_directories(wonton SYSTEM INTERFACE ${FleCSI_INCLUDE_DIR}) + target_include_directories(wonton SYSTEM INTERFACE ${FleCSISP_INCLUDE_DIR}) target_link_libraries(wonton INTERFACE ${FleCSI_LIBRARIES}) target_link_libraries(wonton INTERFACE ${FleCSISP_LIBRARIES}) @@ -341,12 +341,6 @@ if (WONTON_ENABLE_LAPACKE) enable_language(Fortran) include(FortranCInterface) # will ensure the fortran library is linked in - - target_include_directories(wonton INTERFACE ${LAPACKE_INCLUDE_DIRS}) - target_compile_definitions(wonton INTERFACE WONTON_HAS_LAPACKE) - - target_link_libraries(wonton INTERFACE ${LAPACKE_LIBRARIES}) - message(STATUS "LAPACKE_FOUND ${LAPACKE_FOUND}") message(STATUS "LAPACKE_LIBRARIES ${LAPACKE_LIBRARIES}") diff --git a/README.md b/README.md index d3a733e..f0ee3ba 100755 --- a/README.md +++ b/README.md @@ -148,3 +148,32 @@ cmake \ make -j2 ctest --output-on-failure ``` + +## Snow + +Execute the following from the Wonton root directory: + +```c++ +# machine=sn-fey +module load intel/18.0.5 openmpi/2.1.2 cmake/3.14.0 +NGC=/usr/projects/ngc +NGC_INCLUDE_DIR=$NGC/private/include +JALI_INSTALL_PREFIX=$NGC/private/jali/1.1.4-intel-18.0.5-openmpi-2.1.2 +LAPACKE_INSTALL_PREFIX=$NGC/private/lapack/3.8.0-patched-intel-18.0.5 +mkdir build +cd build +cmake \ + -D CMAKE_PREFIX_PATH=$NGC_INCLUDE_DIR \ + -D CMAKE_BUILD_TYPE=Release \ + -D CMAKE_CXX_FLAGS="-Wall -Werror" \ + -D ENABLE_UNIT_TESTS=True \ + -D WONTON_ENABLE_MPI=True \ + -D ENABLE_JENKINS_OUTPUT=True \ + -D WONTON_ENABLE_Jali=True \ + -D Jali_ROOT:FILEPATH=$JALI_INSTALL_PREFIX \ + -D WONTON_ENABLE_LAPACKE=True \ + -D LAPACKE_ROOT:FILEPATH=$LAPACKE_INSTALL_PREFIX \ + .. +make -j36 +ctest -j36 --output-on-failure +``` diff --git a/googletest b/googletest index 549c5d0..2fe3bd9 160000 --- a/googletest +++ b/googletest @@ -1 +1 @@ -Subproject commit 549c5d061e3c5ea4acffb43c92b048dd16812805 +Subproject commit 2fe3bd994b3189899d93f1d5a881e725e046fdc2 diff --git a/jenkins/axis_hpc.yml b/jenkins/axis_hpc.yml new file mode 100644 index 0000000..45113be --- /dev/null +++ b/jenkins/axis_hpc.yml @@ -0,0 +1,18 @@ +COMPILER: + - intel18 + - gcc6 + - gcc7 + +BUILD_TYPE: + - base + - debug + - serial + - readme + - thrust + +exclude: + - COMPILER: gcc6 + BUILD_TYPE: readme + - COMPILER: gcc7 + BUILD_TYPE: readme + diff --git a/jenkins/axis_hpc_pr.yml b/jenkins/axis_hpc_pr.yml new file mode 100644 index 0000000..8d692a8 --- /dev/null +++ b/jenkins/axis_hpc_pr.yml @@ -0,0 +1,15 @@ +COMPILER: + - intel18 + - gcc7 + +BUILD_TYPE: + - base + - debug + - serial + - readme + - thrust + +exclude: + - COMPILER: gcc7 + BUILD_TYPE: readme + diff --git a/jenkins/build_matrix_entry.sh b/jenkins/build_matrix_entry.sh index 6a0be0a..b27a5cd 100755 --- a/jenkins/build_matrix_entry.sh +++ b/jenkins/build_matrix_entry.sh @@ -16,11 +16,27 @@ compiler=$1 build_type=$2 echo "inside build_matrix entry" + # special case for README builds if [[ $build_type == "readme" ]]; then python2 $WORKSPACE/jenkins/parseREADME.py $WORKSPACE/README.md $WORKSPACE exit fi +# special case for README builds +if [[ $build_type == "readme" ]]; then + + # Put a couple of settings in place to generate test output even if + # the README doesn't ask for it. + export CTEST_OUTPUT_ON_FAILURE=1 + CACHE_OPTIONS="-D ENABLE_JENKINS_OUTPUT=True" + sed "s/^ *cmake/& $CACHE_OPTIONS/g" $WORKSPACE/README.md >$WORKSPACE/README.md.1 + python2 $WORKSPACE/jenkins/parseREADME.py \ + $WORKSPACE/README.md.1 \ + $WORKSPACE \ + varan + exit + +fi # set modules and install paths @@ -35,33 +51,31 @@ thrust_dir=${ngc_include_dir} # compiler-specific settings if [[ $compiler == "intel18" ]]; then - intel_version=18.0.1 - cxxmodule=intel/${intel_version} - # openmpi version that libs were built against - openmpi_version=2.1.2 - # openmpi module for compiling and linking - mpi_module=openmpi/2.1.2 - jali_install_dir=$NGC/private/jali/${jali_version}-intel-${intel_version}-openmpi-${openmpi_version} - lapacke_dir=$NGC/private/lapack/${lapack_version}-patched-intel-${intel_version} -elif [[ $compiler == "gcc6" ]]; then - gcc_version=6.4.0 - cxxmodule=gcc/${gcc_version} - # openmpi version that libs were built against + compiler_version=18.0.1 + cxxmodule=intel/${compiler_version} + compiler_suffix="-intel-${compiler_version}" + openmpi_version=2.1.2 - # openmpi module for compiling and linking mpi_module=openmpi/2.1.2 - jali_install_dir=$NGC/private/jali/${jali_version}-gcc-${gcc_version}-openmpi-${openmpi_version} - lapacke_dir=$NGC/private/lapack/${lapack_version}-patched-gcc-${gcc_version} -elif [[ $compiler == "gcc7" ]]; then - gcc_version=7.3.0 - cxxmodule=gcc/${gcc_version} - # openmpi version that libs were built against + mpi_suffix="-openmpi-${openmpi_version}" + +elif [[ $compiler =~ "gcc" ]]; then + + if [[ $compiler == "gcc6" ]]; then + compiler_version=6.4.0 + elif [[ $compiler == "gcc7" ]]; then + compiler_version=7.3.0 + fi + cxxmodule=gcc/${compiler_version} + compiler_suffix="-gcc-${compiler_version}" + openmpi_version=2.1.2 - # openmpi module for compiling and linking - mpi_module=openmpi/2.1.2 - jali_install_dir=$NGC/private/jali/${jali_version}-gcc-${gcc_version}-openmpi-${openmpi_version} - lapacke_dir=$NGC/private/lapack/${lapack_version}-patched-gcc-${gcc_version} + mpi_module=openmpi/${openmpi_version} + mpi_suffix="-openmpi-${openmpi_version}" + fi +jali_install_dir=$NGC/private/jali/${jali_version}${compiler_suffix}${mpi_suffix} +lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} # build-type-specific settings @@ -80,6 +94,13 @@ elif [[ $build_type == "thrust" ]]; then thrust_flags="-D WONTON_ENABLE_THRUST=True -DTHRUST_ROOT:FILEPATH=${thrust_dir}" fi +flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=False" +if [[ $compiler == "gcc6" && $build_type != "serial" ]]; then + flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 + flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 + flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" +fi + export SHELL=/bin/sh export MODULEPATH="" @@ -105,8 +126,9 @@ cmake \ -D ENABLE_JENKINS_OUTPUT=True \ $mpi_flags \ $jali_flags \ + $flecsi_flags \ $lapacke_flags \ $thrust_flags \ .. make -j2 -ctest --output-on-failure +ctest -j2 --output-on-failure diff --git a/jenkins/build_matrix_entry_hpc.sh b/jenkins/build_matrix_entry_hpc.sh new file mode 100755 index 0000000..acd8725 --- /dev/null +++ b/jenkins/build_matrix_entry_hpc.sh @@ -0,0 +1,127 @@ +#!/usr/bin/env bash +# This script is executed on Jenkins using +# +# $WORKSPACE/jenkins/build_matrix_entry.sh +# +# The exit code determines if the test succeeded or failed. +# Note that the environment variable WORKSPACE must be set (Jenkins +# will do this automatically). + +# Exit on error +set -e +# Echo each command +set -x + +compiler=$1 +build_type=$2 + +echo "inside build_matrix entry" + +# special case for README builds +if [[ $build_type == "readme" ]]; then + + # Put a couple of settings in place to generate test output even if + # the README doesn't ask for it. + export CTEST_OUTPUT_ON_FAILURE=1 + CACHE_OPTIONS="-D ENABLE_JENKINS_OUTPUT=True" + sed "s/^ *cmake/& $CACHE_OPTIONS/g" $WORKSPACE/README.md >$WORKSPACE/README.md.1 + python2 $WORKSPACE/jenkins/parseREADME.py \ + $WORKSPACE/README.md.1 \ + $WORKSPACE \ + sn-fey + exit + +fi + +# set modules and install paths + +jali_version=1.1.4 +lapack_version=3.8.0 + +export NGC=/usr/projects/ngc +ngc_include_dir=$NGC/private/include + +thrust_dir=${ngc_include_dir} + + +# compiler-specific settings +if [[ $compiler == "intel18" ]]; then + compiler_version=18.0.5 + cxxmodule=intel/${compiler_version} + compiler_suffix="-intel-${compiler_version}" + + openmpi_version=2.1.2 + mpi_module=openmpi/2.1.2 + mpi_suffix="-openmpi-${openmpi_version}" + +elif [[ $compiler =~ "gcc" ]]; then + + if [[ $compiler == "gcc6" ]]; then + compiler_version=6.4.0 + elif [[ $compiler == "gcc7" ]]; then + compiler_version=7.4.0 + fi + cxxmodule=gcc/${compiler_version} + compiler_suffix="-gcc-${compiler_version}" + + openmpi_version=2.1.2 + mpi_module=openmpi/${openmpi_version} + mpi_suffix="-openmpi-${openmpi_version}" + +fi +jali_install_dir=$NGC/private/jali/${jali_version}${compiler_suffix}${mpi_suffix} +lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} + + +# build-type-specific settings +cmake_build_type=Release +mpi_flags="-D WONTON_ENABLE_MPI=True" +extra_flags= +thrust_flags= +jali_flags="-D WONTON_ENABLE_Jali=True -D Jali_ROOT:FILEPATH=$jali_install_dir" +lapacke_flags="-D WONTON_ENABLE_LAPACKE=True -D LAPACKE_ROOT:FILEPATH=$lapacke_dir" +if [[ $build_type == "debug" ]]; then + cmake_build_type=Debug +elif [[ $build_type == "serial" ]]; then + mpi_flags= + jali_flags= # jali is not available in serial +elif [[ $build_type == "thrust" ]]; then + thrust_flags="-D WONTON_ENABLE_THRUST=True -DTHRUST_ROOT:FILEPATH=${thrust_dir}" +fi + +# We don't seem to have a suitable (read as "ancient") flecsi install on HPC +#if [[ $compiler == "gcc6" && $build_type != "serial" ]]; then +# flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 +# flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 +# flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" +#fi + +export SHELL=/bin/sh +#Rely on default user environment to load modules; these scripts can be found in /etc/profile.d/ as of 8/25/20 the scripts that set up modules on snow are /etc/profile.d/z00_lmod.sh; /etc/profile.d/00-modulepath.sh; /etc/profile.d/z01-modules.lanl.sh; +module load $cxxmodule +if [[ -n "$mpi_flags" ]]; then + module load ${mpi_module} +fi +module load cmake/3.14.0 # 3.13 or higher is required + +echo $WORKSPACE +cd $WORKSPACE + +mkdir build +cd build + +cmake \ + -D CMAKE_BUILD_TYPE=$cmake_build_type \ + -D CMAKE_CXX_FLAGS="-Wall -Werror" \ + -D CMAKE_PREFIX_PATH=$ngc_include_dir \ + -D ENABLE_UNIT_TESTS=True \ + -D ENABLE_APP_TESTS=True \ + -D ENABLE_JENKINS_OUTPUT=True \ + $mpi_flags \ + $jali_flags \ + $flecsi_flags \ + $lapacke_flags \ + $thrust_flags \ + .. +make -j36 +ctest -j36 --output-on-failure diff --git a/jenkins/install_matrix_entry.sh b/jenkins/install_matrix_entry.sh index 6248744..188951e 100755 --- a/jenkins/install_matrix_entry.sh +++ b/jenkins/install_matrix_entry.sh @@ -26,7 +26,7 @@ export NGC=/usr/local/codes/ngc ngc_include_dir=$NGC/private/include wonton_install_dir=$NGC/private/wonton -wonton_version=dev +wonton_version=$3 thrust_dir=${ngc_include_dir} @@ -67,9 +67,11 @@ jali_flags="-D WONTON_ENABLE_Jali:BOOL=True -D Jali_ROOT:PATH=$jali_install_dir" lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} lapacke_flags="-D WONTON_ENABLE_LAPACKE:BOOL=True -D LAPACKE_ROOT:PATH=$lapacke_dir" +flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=False" if [[ $compiler == "gcc6" ]]; then flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 - flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir" + flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 + flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" fi diff --git a/jenkins/parseREADME.py b/jenkins/parseREADME.py index 81eea5d..6748b26 100755 --- a/jenkins/parseREADME.py +++ b/jenkins/parseREADME.py @@ -9,11 +9,16 @@ # line of a block be a followed by a hostname. The script will then # execute all commands specified after that point on the host machine # until the block is closed with a ```. -# Args: Arg 1 is file to parse, Arg 2 is the workspace +# Args: Arg 1 is file to parse, Arg 2 is the workspace, Arg 3 is the machine + +infile = sys.argv[1] +workspace = sys.argv[2] +machine = sys.argv[3] + code = False scripts = {} mach_kw = "machine=" -for line in open(sys.argv[1]): +for line in open(infile): if line.startswith("```"): code = not code currentScript = None @@ -24,21 +29,19 @@ # with the "machine=" keyword if not line.startswith("#") and mach_kw not in line: continue host = line.split(mach_kw)[-1].strip() - currentScript = ["set -e\n pwd \n cd "+ sys.argv[2] + "\n rm -Rf build \n"] + currentScript = ["set -e\n pwd \n cd "+ workspace + "\n rm -Rf build \n"] scripts[host] = currentScript else: currentScript.append(line) retcode = 0 for host in scripts: - print(host) - if(host.startswith('varan')): + if(host == machine): print "HOST " + host - scripts[host] += 'cd ..;rm -Rf build;' print scripts[host] sshHost = host.split(':')[0] proc = subprocess.Popen(['bash'], stdin = subprocess.PIPE) print "opened process" proc.communicate("".join(scripts[host])) if proc.returncode: - retcode = proc.returncode + retcode = proc.returncode sys.exit(retcode) diff --git a/wonton/mesh/AuxMeshTopology.h b/wonton/mesh/AuxMeshTopology.h index e207e3f..54e588f 100644 --- a/wonton/mesh/AuxMeshTopology.h +++ b/wonton/mesh/AuxMeshTopology.h @@ -18,7 +18,9 @@ Please see the license file at the root of this repository, or at: #include #include "wonton/support/wonton.h" +#include "wonton/support/CoordinateSystem.h" #include "wonton/support/Point.h" +#include "wonton/support/Polytope.h" namespace Wonton { @@ -505,7 +507,7 @@ class AuxMeshTopology { template void cell_centroid(int const cellid, Point *ccen) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif for (int d = 0; d < D; ++d) @@ -515,7 +517,7 @@ class AuxMeshTopology { //! Volume of a cell double cell_volume(int const cellid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif assert(sides_requested_); @@ -527,7 +529,7 @@ class AuxMeshTopology { template void face_centroid(int const faceid, Point *fcen) const { -#ifdef DEBUG +#ifndef NDEBUG assert(faceid < num_entities(FACE, ALL)); #endif for (int d = 0; d < D; ++d) @@ -545,7 +547,7 @@ class AuxMeshTopology { //! centroid form a positive volume tet. int side_get_node(int const sideid, int const inode) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(sides_requested_); @@ -557,7 +559,7 @@ class AuxMeshTopology { //! Cell of side int side_get_cell(int const sideid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(sides_requested_); @@ -568,7 +570,7 @@ class AuxMeshTopology { //! Face of side int side_get_face(int const sideid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(sides_requested_); @@ -584,7 +586,7 @@ class AuxMeshTopology { //! So, side_get_node(s,i) = wedge_get_node(side_get_wedge(s,i)) int side_get_wedge(int const sideid, int iwedge) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(wedges_requested_); @@ -600,7 +602,7 @@ class AuxMeshTopology { //! returns -1 int side_get_opposite_side(int const sideid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(sides_requested_); @@ -611,7 +613,7 @@ class AuxMeshTopology { //! Get all the sides of a cell void cell_get_sides(int const cellid, std::vector *csides) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif assert(sides_requested_); @@ -667,7 +669,7 @@ class AuxMeshTopology { //! Volume of a side double side_volume(int const sideid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(sideid < num_entities(SIDE, ALL)); #endif assert(sides_requested_); @@ -677,7 +679,7 @@ class AuxMeshTopology { //! Side of wedge int wedge_get_side(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(wedges_requested_); @@ -688,7 +690,7 @@ class AuxMeshTopology { //! Cell of wedge int wedge_get_cell(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(wedges_requested_); @@ -700,7 +702,7 @@ class AuxMeshTopology { //! Face of wedge int wedge_get_face(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(wedges_requested_); @@ -712,7 +714,7 @@ class AuxMeshTopology { //! Corner of a wedge int wedge_get_corner(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(corners_requested_); @@ -723,7 +725,7 @@ class AuxMeshTopology { //! node of a wedge int wedge_get_node(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(wedges_requested_); @@ -740,7 +742,7 @@ class AuxMeshTopology { //! routine returns -1 int wedge_get_opposite_wedge(const int wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif assert(wedges_requested_); @@ -760,7 +762,7 @@ class AuxMeshTopology { //! 2D int wedge_get_adjacent_wedge(const int wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif // Wedges come in pairs; their IDs are (2*sideid) and (2*sideid+1) @@ -775,7 +777,7 @@ class AuxMeshTopology { //! Volume of a wedge - half its side volume double wedge_volume(int const wedgeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(wedgeid < num_entities(WEDGE, ALL)); #endif int sideid = static_cast(wedgeid/2); @@ -832,7 +834,7 @@ class AuxMeshTopology { //! Get all the wedges in a cell void cell_get_wedges(int const cellid, std::vector *wedgeids) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif int nsides = cell_side_ids_[cellid].size(); @@ -849,7 +851,7 @@ class AuxMeshTopology { void node_get_wedges(int const nodeid, Entity_type const type, std::vector *wedgeids) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(wedges_requested_ && corners_requested_); @@ -871,7 +873,7 @@ class AuxMeshTopology { //! Get node of corner int corner_get_node(const int cornerid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cornerid < num_entities(CORNER, ALL)); #endif assert(corners_requested_); @@ -882,7 +884,7 @@ class AuxMeshTopology { //! Get cell of corner int corner_get_cell(int const cornerid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cornerid < num_entities(CORNER, ALL)); #endif assert(corners_requested_); @@ -893,7 +895,7 @@ class AuxMeshTopology { //! Get wedges of a corner void corner_get_wedges(int const cornerid, std::vector *wedgeids) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cornerid < num_entities(CORNER, ALL)); #endif assert(corners_requested_); @@ -907,7 +909,7 @@ class AuxMeshTopology { void node_get_corners(int const nodeid, Entity_type const type, std::vector *cornerids) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(corners_requested_); @@ -930,7 +932,7 @@ class AuxMeshTopology { //! Get corners in a cell void cell_get_corners(int const cellid, std::vector *cornerids) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif assert(corners_requested_); @@ -942,7 +944,7 @@ class AuxMeshTopology { //! Get a cell's corner at a particular node of the cell int cell_get_corner_at_node(int const cellid, int const nodeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif assert(corners_requested_); @@ -957,7 +959,7 @@ class AuxMeshTopology { //! Volume of a corner double corner_volume(int const cornerid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cornerid < num_entities(CORNER, ALL)); #endif assert(corners_requested_); @@ -985,7 +987,7 @@ class AuxMeshTopology { void decompose_cell_into_tets(int cellid, std::vector, 4>> *tcoords, const bool planar_hex) const { -#ifdef DEBUG +#ifndef NDEBUG assert(cellid < num_entities(CELL, ALL)); #endif if (planar_hex @@ -1051,7 +1053,7 @@ class AuxMeshTopology { void dual_cell_get_coordinates(int const nodeid, std::vector> *pplist) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(basicmesh_ptr_->space_dimension() == 2); @@ -1186,7 +1188,7 @@ class AuxMeshTopology { void dual_cell_get_node_adj_cells(int const nodeid, Entity_type const ptype, std::vector *adjnodes) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif basicmesh_ptr_->node_get_cell_adj_nodes(nodeid, ptype, adjnodes); @@ -1200,7 +1202,7 @@ class AuxMeshTopology { void dual_cell_get_coordinates(int const nodeid, std::vector> *pplist) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(basicmesh_ptr_->space_dimension() == 3); @@ -1262,7 +1264,7 @@ class AuxMeshTopology { void dual_wedges_get_coordinates(int nodeid, std::vector, 4>> *wcoords) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(wedges_requested_); @@ -1281,7 +1283,7 @@ class AuxMeshTopology { template void dual_cell_centroid(int nodeid, Point *centroid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(wedges_requested_); @@ -1314,7 +1316,7 @@ class AuxMeshTopology { //! Get the volume of dual cell by finding the corners that attach to the node double dual_cell_volume(int const nodeid) const { -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(corners_requested_); @@ -1328,6 +1330,7 @@ class AuxMeshTopology { protected: + template void build_aux_entities() { int dim = basicmesh_ptr_->space_dimension(); @@ -1343,58 +1346,52 @@ class AuxMeshTopology { // // So the steps should be: // + // compute_cell_moments // compute_approximate_face_centroids - // compute_approximate_cell_centroids // // build_sides (topology only) // // compute_face_centroids - // compute_cell_centroids // // compute_side_volumes if (dim == 1) { + compute_cell_moments<1, CoordSys>(); compute_approximate_face_centroids<1>(); - compute_approximate_cell_centroids<1>(); if (sides_requested_) build_sides_1D(*this); compute_face_centroids<1>(); - compute_cell_centroids<1>(); if (sides_requested_) - compute_side_volumes<1>(); + compute_side_volumes<1, CoordSys>(); } else if (dim == 2) { + compute_cell_moments<2, CoordSys>(); compute_approximate_face_centroids<2>(); - compute_approximate_cell_centroids<2>(); if (sides_requested_) build_sides_2D(*this); compute_face_centroids<2>(); - compute_cell_centroids<2>(); if (sides_requested_) - compute_side_volumes<2>(); + compute_side_volumes<2, CoordSys>(); } else if (dim == 3) { + compute_cell_moments<3, CoordSys>(); compute_approximate_face_centroids<3>(); - compute_approximate_cell_centroids<3>(); if (sides_requested_) build_sides_3D(*this); compute_face_centroids<3>(); - compute_cell_centroids<3>(); if (sides_requested_) - compute_side_volumes<3>(); + compute_side_volumes<3, CoordSys>(); } - compute_cell_volumes(); // needs side volumes - if (wedges_requested_) build_wedges(); if (corners_requested_) build_corners(); @@ -1406,10 +1403,10 @@ class AuxMeshTopology { private: template void compute_approximate_cell_centroids(); template void compute_approximate_face_centroids(); - template void compute_cell_centroids(); - template void compute_face_centroids(); - template void compute_side_volumes(); - void compute_cell_volumes(); + template void compute_face_centroids(); // should be templated later on CoordSys + + template void compute_cell_moments(); + template void compute_side_volumes(); void build_wedges(); void build_corners(); @@ -2224,6 +2221,7 @@ void AuxMeshTopology::compute_face_centroids() { // Compute an approximate centroid as the geometric mean of the cell nodes +// NOT USED - CAN BE REMOVED IN THE NEAR FUTURE template template @@ -2252,58 +2250,103 @@ void AuxMeshTopology::compute_approximate_cell_centroids() { } } -// Use approximate cell centroids and the geometry of sides (tets) -// using these approximate centroids to compute the real -// centroids. The real centroid is computed as the volume weighted -// average of the centroids of the sides (tets) + +// Compute centroids and volumes for cells template -template -void AuxMeshTopology::compute_cell_centroids() { +template +void AuxMeshTopology::compute_cell_moments() { int ncells = basicmesh_ptr_->num_owned_cells() + basicmesh_ptr_->num_ghost_cells(); - // Resize (just to be sure) but don't reinitialize since it contains - // the approximate centroids (geometric centers) + cell_volumes_.clear(); + cell_volumes_.assign(ncells, 0.0); - cell_centroids_.resize(ncells); + cell_centroids_.clear(); + std::vector pnt(dim, 0.0); + cell_centroids_.assign(ncells, pnt); std::array, dim+1> sxyz; - bool posvol_order = true; for (int c = 0; c < ncells; ++c) { - // Now compute the real centroid as the volume weighted average of the - // centroids of the sides with the temporary geometry - - std::vector csides; - basicmesh_ptr_->cell_get_sides(c, &csides); + // fork the code since construction of the Polytope + // object in 1D and 3D requires more discussion - Point ccen; double cellvol = 0.0; - for (auto const& s : csides) { - side_get_coordinates(s, &sxyz, posvol_order); - Point scen; - for (int i = 0; i < dim+1; i++) - scen += sxyz[i]; - scen /= (dim+1); + Point ccen; + + std::vector cnodes; + basicmesh_ptr_->cell_get_nodes(c, &cnodes); + int ncnodes = cnodes.size(); - double svol = calc_side_volume(sxyz); - cellvol += svol; + // build list of face nodes (not really needed for 2D and 1D but + // we have to do it for consistency of the Polytope constructor + // further below). + + std::vector cfaces, cfdirs; + basicmesh_ptr_->cell_get_faces_and_dirs(c, &cfaces, &cfdirs); + int nfaces = cfaces.size(); + + // map of face nodes into cell node list + std::vector> fnmap(nfaces); + + for (int i = 0; i < nfaces; i++) { + std::vector fnodes; + basicmesh_ptr_->face_get_nodes(cfaces[i], &fnodes); - ccen += scen*svol; + int nfnodes = fnodes.size(); + int dir = cfdirs[i]; + + std::vector fnodes1(nfnodes); + for (int j = 0; j < nfnodes; j++) { + int n = (dir == 1) ? fnodes[j] : fnodes[nfnodes-j-1]; + for (int k = 0; k < ncnodes; k++) + if (cnodes[k] == n) { + fnodes1[j] = k; + break; + } + } + fnmap[i] = fnodes1; } + + std::vector> cpoints(cnodes.size()); + for (int i = 0; i < ncnodes; i++) + basicmesh_ptr_->node_get_coordinates(cnodes[i], &(cpoints[i])); + + Polytope poly(cpoints, fnmap); + + int order = 1; + if (std::is_same::value) + order = 2; - //If we have a degenerate cell of zero volume, we use the geometric center + auto moments = poly.moments(order); + CoordSys::template shift_moments_list(moments); + + cellvol = moments[0]; + for (int d = 0; d < dim; d++) + ccen[d] = moments[d + 1]; + if (cellvol > std::numeric_limits::epsilon()) { ccen /= cellvol; for (int d = 0; d < dim; d++) cell_centroids_[c][d] = ccen[d]; // true centroid + } else { + // use the geometric center of cell nodes as a proxy for the centroid + Point geomcen; + for (int n = 0; n < ncnodes; ++n) { + geomcen += cpoints[n]; + } + geomcen /= ncnodes; + for (int d = 0; d < dim; ++d) + cell_centroids_[c][d] = geomcen[d]; } + + cell_volumes_[c] = cellvol; } } template -template +template void AuxMeshTopology::compute_side_volumes() { int num_sides_all = num_sides_owned_ + num_sides_ghost_; side_volumes_.resize(num_sides_all); @@ -2313,23 +2356,16 @@ void AuxMeshTopology::compute_side_volumes() { for (int s = 0; s < num_sides_all; s++) { side_get_coordinates(s, &sxyz, posvol_order); side_volumes_[s] = calc_side_volume(sxyz); + + // sufficient for correct cell-volume calculation in curvilinear coordinate systems + // volume corresponds to 1 radiant of the cylindrically symmetric side + if (std::is_same::value) { + side_volumes_[s] *= (sxyz[0][0] + sxyz[1][0] + sxyz[2][0]) / 3; + } } } -template -void AuxMeshTopology::compute_cell_volumes() { - int ncells = basicmesh_ptr_->num_owned_cells() + - basicmesh_ptr_->num_ghost_cells(); - - cell_volumes_.clear(); - cell_volumes_.assign(ncells, 0.0); - - for (int c = 0; c < ncells; ++c) - for (auto s : cell_side_ids_[c]) - cell_volumes_[c] += side_volumes_[s]; -} - //! coords of nodes of a cell template template @@ -2443,7 +2479,7 @@ void AuxMeshTopology::dual_cell_get_facetization(int const nodeid, facetpoints->clear(); points->clear(); -#ifdef DEBUG +#ifndef NDEBUG assert(nodeid < num_entities(NODE, ALL)); #endif assert(wedges_requested_); diff --git a/wonton/mesh/direct_product/direct_product_mesh.h b/wonton/mesh/direct_product/direct_product_mesh.h index aa0ada8..b757d5c 100644 --- a/wonton/mesh/direct_product/direct_product_mesh.h +++ b/wonton/mesh/direct_product/direct_product_mesh.h @@ -448,7 +448,7 @@ double Direct_Product_Mesh::get_axis_point( assert(dim >= 0); assert(dim < D); -#if DEBUG +#ifndef NDEBUG int npall = axis_points_[dim].size(); assert(pointid >= -num_ghost_layers_); assert(pointid <= npall); diff --git a/wonton/mesh/flecsi/flecsi_mesh_wrapper.h b/wonton/mesh/flecsi/flecsi_mesh_wrapper.h index d748db1..94fbbb8 100644 --- a/wonton/mesh/flecsi/flecsi_mesh_wrapper.h +++ b/wonton/mesh/flecsi/flecsi_mesh_wrapper.h @@ -119,11 +119,13 @@ class flecsi_mesh_t : public AuxMeshTopology> { //! \param[in] mesh The minimum coordinates of the domain. explicit flecsi_mesh_t(const mesh_t & mesh, bool request_sides = true, - bool request_wedges = true, - bool request_corners = true) : - cells_(mesh.cells()), faces_(mesh.faces()), - vertices_(mesh.vertices()), mesh_(&mesh), - wonton_mesh_aux_t(request_sides, request_wedges, request_corners) { + bool request_wedges = true, + bool request_corners = true) + : wonton_mesh_aux_t(request_sides, request_wedges, request_corners), + mesh_(&mesh), + cells_(mesh.cells()), + faces_(mesh.faces()), + vertices_(mesh.vertices()) { // base class (AuxMeshTopology) method that has to be called here // and not in the constructor of the base class because it needs // access to methods in this class which in turn need access to diff --git a/wonton/mesh/flecsi/test/flecsi-headers.h b/wonton/mesh/flecsi/test/flecsi-headers.h new file mode 100644 index 0000000..12d13fb --- /dev/null +++ b/wonton/mesh/flecsi/test/flecsi-headers.h @@ -0,0 +1,16 @@ +/* + * This file is part of the Ristra Wonton project. + * Please see the license file at the root of this repository, or at: + * https://github.com/laristra/wonton/blob/master/LICENSE + */ + +#ifndef WONTON_FLECSI_HEADERS_H +#define WONTON_FLECSI_HEADERS_H + +#pragma GCC system_header + +#include "flecsi-sp.h" +#include "flecsi-sp/burton/burton.h" +#include "flecsi-sp/burton/factory.h" + +#endif //WONTON_FLECSI_HEADERS_H diff --git a/wonton/mesh/flecsi/test/test_flecsi_mesh_wrapper.cc b/wonton/mesh/flecsi/test/test_flecsi_mesh_wrapper.cc index 068ce4e..636d4e8 100644 --- a/wonton/mesh/flecsi/test/test_flecsi_mesh_wrapper.cc +++ b/wonton/mesh/flecsi/test/test_flecsi_mesh_wrapper.cc @@ -11,10 +11,7 @@ Please see the license file at the root of this repository, or at: #include "gtest/gtest.h" #include "mpi.h" -#include "flecsi-sp.h" -#include "flecsi-sp/burton/burton.h" -#include "flecsi-sp/burton/factory.h" - +#include "wonton/mesh/flecsi/test/flecsi-headers.h" #include "wonton/support/Point.h" namespace fmesh = flecsi::sp::burton; @@ -603,7 +600,7 @@ TEST(FleCSI_Mesh_Wrapper, Decompose_Cell_Into_Tets) { TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_2D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc > 1) return; #endif @@ -620,8 +617,6 @@ TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_2D) { bool request_corners = true; flecsi_mesh_2d_t mesh_wrapper(mesh, request_sides, request_wedges, request_corners); - double dp; - int ncells = mesh_wrapper.num_entities(Wonton::CELL, Wonton::ALL); for (int c = 0; c < ncells; ++c) { std::vector csides; @@ -738,7 +733,7 @@ TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_2D) { TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_3D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc > 1) return; #endif @@ -759,7 +754,6 @@ TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_3D) { flecsi_mesh_3d_t mesh_wrapper(mesh, request_sides, request_wedges, request_corners); int ncells = mesh_wrapper.num_entities(Wonton::CELL, Wonton::ALL); - double dp; for (int c = 0; c < ncells; ++c) { std::vector csides; mesh_wrapper.cell_get_sides(c, &csides); @@ -881,7 +875,7 @@ TEST(FleCSI_Mesh_Wrapper, MESH_SIDES_3D) { TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_2D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc > 1) return; #endif @@ -899,7 +893,6 @@ TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_2D) { bool request_corners = true; flecsi_mesh_2d_t mesh_wrapper(mesh, request_sides, request_wedges, request_corners); - double dp; int ncells = mesh_wrapper.num_entities(Wonton::CELL, Wonton::ALL); for (int c = 0; c < ncells; ++c) { std::vector cwedges; @@ -990,8 +983,6 @@ TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_2D) { std::array, 3> wcoords2; mesh_wrapper.wedge_get_coordinates(w, &wcoords2, true); - bool flipped = true; - if (wcoords[1][0] == wcoords2[2][0] && wcoords[1][1] == wcoords2[2][1] && wcoords[2][0] == wcoords2[1][0] && @@ -1028,7 +1019,7 @@ TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_2D) { TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_3D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc > 1) return; #endif @@ -1049,7 +1040,6 @@ TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_3D) { flecsi_mesh_3d_t mesh_wrapper(mesh, request_sides, request_wedges, request_corners); - double dp; int ncells = mesh_wrapper.num_entities(Wonton::CELL, Wonton::ALL); for (int c = 0; c < ncells; ++c) { std::vector cwedges; @@ -1182,7 +1172,7 @@ TEST(FleCSI_Mesh_Wrapper, MESH_WEDGES_3D) { TEST(FleCSI_Mesh_Wrapper, MESH_CORNERS_2D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc > 1 ) return; #endif @@ -1302,7 +1292,7 @@ TEST(FleCSI_Mesh_Wrapper, MESH_CORNERS_2D) { TEST(FleCSI_Mesh_Wrapper, MESH_CORNERS_3D) { #ifdef WONTON_ENABLE_MPI - int nproc, me; + int nproc; MPI_Comm_size(MPI_COMM_WORLD, &nproc); if (nproc>1) return; #endif diff --git a/wonton/mesh/simple/simple_mesh_wrapper.h b/wonton/mesh/simple/simple_mesh_wrapper.h index ab0e6e2..7586757 100644 --- a/wonton/mesh/simple/simple_mesh_wrapper.h +++ b/wonton/mesh/simple/simple_mesh_wrapper.h @@ -13,6 +13,7 @@ Please see the license file at the root of this repository, or at: #include #include "wonton/mesh/AuxMeshTopology.h" +#include "wonton/support/CoordinateSystem.h" #include "wonton/support/wonton.h" #include "wonton/support/Point.h" @@ -47,14 +48,22 @@ class Simple_Mesh_Wrapper : public AuxMeshTopology { @param[in] request_corners Should the AuxMeshToplogy class build corner datastructures? */ + + // Typically, we template a class to add functionality. By majority vote, + // this has not been done for curvilinier coordinate systems. Instead, + // a low-intrusion (optional parameter) approach was adopted explicit Simple_Mesh_Wrapper(Simple_Mesh const & mesh, bool request_sides = true, bool request_wedges = true, - bool request_corners = true) : + bool request_corners = true, + CoordSysType coord_sys = CoordSysType::Cartesian) : AuxMeshTopology( request_sides, request_wedges, request_corners), mesh_(mesh) { - AuxMeshTopology::build_aux_entities(); + if (coord_sys == CoordSysType::CylindricalAxisymmetric) + AuxMeshTopology::build_aux_entities(); + else + AuxMeshTopology::build_aux_entities(); } // explicit constructor /// Copy constructor (disabled). diff --git a/wonton/mesh/simple/test/test_simple_mesh_wrapper.cc b/wonton/mesh/simple/test/test_simple_mesh_wrapper.cc index 88563cd..37b0d10 100644 --- a/wonton/mesh/simple/test/test_simple_mesh_wrapper.cc +++ b/wonton/mesh/simple/test/test_simple_mesh_wrapper.cc @@ -11,6 +11,7 @@ Please see the license file at the root of this repository, or at: #include "wonton/mesh/simple/simple_mesh.h" #include "wonton/mesh/simple/simple_mesh_wrapper.h" +#include "wonton/support/CoordinateSystem.h" #include "wonton/support/wonton.h" #include "wonton/support/Point.h" #include "wonton/support/Vector.h" diff --git a/wonton/state/test/test_state_manager.cc b/wonton/state/test/test_state_manager.cc index 7170f27..903a87d 100644 --- a/wonton/state/test/test_state_manager.cc +++ b/wonton/state/test/test_state_manager.cc @@ -371,9 +371,9 @@ TEST(StateManager, multiMatField){ // test that the value obtained through the state manager is correct for (auto& kv : out->get_data()) { - auto& out=kv.second; - for (unsigned j=0; j < out.size(); j++){ - ASSERT_EQ(data[kv.first][j], out[j]); + auto& tmp=kv.second; + for (unsigned j=0; j < tmp.size(); j++){ + ASSERT_EQ(data[kv.first][j], tmp[j]); } } @@ -404,9 +404,9 @@ TEST(StateManager, mixedFields){ // test that the value obtained through the state manager is correct for (auto& kv : out->get_data()) { - auto& out=kv.second; - for (unsigned j=0; j < out.size(); j++){ - ASSERT_EQ(data[kv.first][j], out[j]); + auto& tmp=kv.second; + for (unsigned j=0; j < tmp.size(); j++){ + ASSERT_EQ(data[kv.first][j], tmp[j]); } } diff --git a/wonton/support/CMakeLists.txt b/wonton/support/CMakeLists.txt index bae6b46..9123597 100644 --- a/wonton/support/CMakeLists.txt +++ b/wonton/support/CMakeLists.txt @@ -43,6 +43,15 @@ if (WONTON_ENABLE_MPI) target_compile_definitions(wonton_support PUBLIC OMPI_SKIP_MPICXX) endif () +#----------------------------------------------------------------------------- +# Lapacke +#----------------------------------------------------------------------------- +if (WONTON_ENABLE_LAPACKE AND LAPACKE_FOUND) + target_include_directories(wonton_support PUBLIC ${LAPACKE_INCLUDE_DIRS}) + target_compile_definitions(wonton_support PUBLIC WONTON_HAS_LAPACKE) + target_link_libraries(wonton_support PUBLIC ${LAPACKE_LIBRARIES}) +endif () + #----------------------------------------------------------------------------- # Thrust information #----------------------------------------------------------------------------- diff --git a/wonton/support/CoordinateSystem.h b/wonton/support/CoordinateSystem.h index c43f86b..ef9f463 100644 --- a/wonton/support/CoordinateSystem.h +++ b/wonton/support/CoordinateSystem.h @@ -74,9 +74,13 @@ namespace CoordinateSystem { // Swap vectors new_moments.swap(moments); } - } +// Typically, we template a class to add functionality. By majority vote, this +// has not been done for curvilinier coordinate systems. +enum class CoordSysType { Cartesian = 1, + CylindricalAxisymmetric }; + // ============================================================================ /// Cartesian Coordinates struct CartesianCoordinates { diff --git a/wonton/support/Point.h b/wonton/support/Point.h index 5b9e024..a87f790 100644 --- a/wonton/support/Point.h +++ b/wonton/support/Point.h @@ -36,10 +36,6 @@ namespace Wonton { -const int X = 0; -const int Y = 1; -const int Z = 2; - /*! @class Point "Point.h" @brief Represents a point in an N-dimensional space. diff --git a/wonton/support/Polytope.h b/wonton/support/Polytope.h index 3a15008..5db5b59 100644 --- a/wonton/support/Polytope.h +++ b/wonton/support/Polytope.h @@ -21,31 +21,41 @@ template class Polytope { public: /*! - @brief Constructor for a 2D polygon + @brief Constructor for a 1D or 2D polygon @param vertex_points Vector of coordinates of the polygon's vertices listed counter-clockwise */ - explicit Polytope(const std::vector< Point<2> >& vertex_points) { + Polytope(const std::vector< Point >& vertex_points) { + assert(D == 1 || D == 2); // this signature won't work for 3D polytopes + int nfaces = static_cast(vertex_points.size()); - assert(nfaces > 2); + assert((D == 1 && nfaces == 2) || (nfaces >= D+1)); vertex_points_ = vertex_points; face_vertices_.resize(nfaces); - for (int iface = 0; iface < nfaces; iface++) - face_vertices_[iface] = { iface, (iface + 1)%nfaces }; + if (D == 1) + for (int iface = 0; iface < nfaces; iface++) + face_vertices_[iface] = { iface }; + else + for (int iface = 0; iface < nfaces; iface++) + face_vertices_[iface] = { iface, (iface + 1)%nfaces }; } /*! - @brief Constructor for a 3D polyhedron + @brief Constructor for a general polyhedron @param vertex_points Vector of coordinates of the polyhedron's vertices - @param face_vertices Faces of the polyhedron, every face is given by - IDs of its vertices in counter-clockwise order (i.e. so that the normal - to the face is pointing outside of the polyhedron) + @param face_vertices Face vertices of the polyhedron + + In 3D, every face is given by IDs of its vertices in + counter-clockwise order (i.e. so that the normal to the face is + pointing outside of the polyhedron). In 2D, a "face" is an edge between two */ - Polytope(const std::vector< Point<3> >& vertex_points, + Polytope(const std::vector< Point >& vertex_points, const std::vector< std::vector >& face_vertices) { - assert(vertex_points.size() > 3); - assert(face_vertices.size() > 3); + assert((D == 1 && vertex_points.size() == 2) || + (vertex_points.size() >= D+1)); + assert((D == 1 && face_vertices.size() == 2) || + (face_vertices.size() >= D+1)); vertex_points_ = vertex_points; face_vertices_ = face_vertices; @@ -123,8 +133,8 @@ class Polytope { /*! @brief Moments for a face of the polytope @param face_id ID of the face of the polytope - @return Vector of moments; moments[0] is the size, moments[i+1]/moments[0] is i-th - coordinate of the centroid + @return Vector of moments; moments[0] is the area (length in 2D), moments[i+1]/moments[0] + is i-th coordinate of the centroid. */ std::vector face_moments(int face_id) const; @@ -137,10 +147,10 @@ class Polytope { /*! @brief Moments for the polytope - @return Vector of moments; moments[0] is the size, moments[i+1]/moments[0] is i-th - coordinate of the centroid + @return Vector of moments; moments[0] is the volume (area in 2D), moments[i+1]/moments[0] + is i-th coordinate of the centroid */ - std::vector moments() const; + std::vector moments(int order = 1) const; private: std::vector< Point > vertex_points_; // coordinates of vertices @@ -286,14 +296,27 @@ Vector<3> Polytope<3>::face_normal(int face_id) const { } /*! - @brief Moments for the polygon + @brief Moments for the polygon (1D) + @param moments Computed moments: moments[0] is length, + moments[i+1]/moments[0] is i-th coordinate of the centroid, i=1 +*/ +template<> +inline +std::vector Polytope<1>::moments(int order) const { + throw std::runtime_error("Wonton::Polytope<1>::moments(int) not implemented"); +} + +/*! + @brief Moments for the polygon (2D) @param moments Computed moments: moments[0] is area, moments[i+1]/moments[0] is i-th coordinate of the centroid, i=1,2 */ template<> inline -std::vector Polytope<2>::moments() const { - std::vector poly_moments(3, 0.0); +std::vector Polytope<2>::moments(int order) const { + assert(order < 3); + int nmoments = (order + 1) * (order + 2) / 2; + std::vector poly_moments(nmoments, 0.0); int nvrts = num_vertices(); for (int ivrt = 0; ivrt < nvrts; ivrt++) { @@ -301,13 +324,24 @@ std::vector Polytope<2>::moments() const { double cur_term = vertex_points_[ifv][0]*vertex_points_[isv][1] - vertex_points_[ifv][1]*vertex_points_[isv][0]; poly_moments[0] += cur_term; - for (int ixy = 0; ixy < 2; ixy++) - poly_moments[ixy + 1] += cur_term*( - vertex_points_[ifv][ixy] + vertex_points_[isv][ixy]); + if (order > 0) { + for (int n = 0; n < 2; ++n) + poly_moments[n + 1] += cur_term * (vertex_points_[ifv][n] + vertex_points_[isv][n]); + } + if (order > 1) { + double xa = vertex_points_[ifv][0], ya = vertex_points_[ifv][1]; + double xb = vertex_points_[isv][0], yb = vertex_points_[isv][1]; + poly_moments[3] += cur_term * (xa * xa + xa * xb + xb * xb); + poly_moments[4] += cur_term * (xa * ya + xb * yb + (xa * yb + xb * ya) / 2); + poly_moments[5] += cur_term * (ya * ya + ya * yb + yb * yb); + } } + poly_moments[0] /= 2.0; - for (int ixy = 0; ixy < 2; ixy++) - poly_moments[ixy + 1] /= 6.0; + for (int n = 0; n < 2; ++n) poly_moments[n + 1] /= 6.0; + if (order > 1) { + for (int n = 0; n < 3; ++n) poly_moments[n + 3] /= 12.0; + } return poly_moments; } @@ -320,7 +354,8 @@ std::vector Polytope<2>::moments() const { */ template<> inline -std::vector Polytope<3>::moments() const { +std::vector Polytope<3>::moments(int order) const { + assert(order < 2); std::vector poly_moments(4, 0.0); int nfaces = num_faces(); diff --git a/wonton/support/equifactor.h b/wonton/support/equifactor.h index 1e579c4..aba9719 100644 --- a/wonton/support/equifactor.h +++ b/wonton/support/equifactor.h @@ -41,17 +41,12 @@ Please see the license file at the root of this repository, or at: namespace Wonton { -#ifdef ENABLE_DEBUG +#ifndef NDEBUG +static // for ODR void print_sets(std::vector> sets); #endif -// gcc 7.3.0 doesn't recognize that this function is used in the code -// so use the compiler specific attribute to turn off the warning (since we -// use -Werror and cannot get past compilation) -#if defined(__GNUC__) && (__GNUC__ == 7 && __GNUC_MINOR__ == 3) -__attribute__ ((unused)) -#endif -static inline +static inline // for ODR std::vector equifactor(int const N, int const D, int const randseed = 0) { clock_t startclock, curclock; startclock = clock(); @@ -85,7 +80,7 @@ std::vector equifactor(int const N, int const D, int const randseed = 0) { kset++; } -#ifdef ENABLE_DEBUG +#ifndef NDEBUG // print out initial sets std::cerr << "\n\nNumber of sets D " << D << "\n"; std::cerr << "INITIAL STATE:\n"; @@ -258,7 +253,7 @@ std::vector equifactor(int const N, int const D, int const randseed = 0) { if (outer_iter > maxiter) outer_done = true; } // while (!outer_done) -#ifdef ENABLE_DEBUG +#ifndef NDEBUG std::cerr << "\nFINAL STATE:\n"; std::vector> newsets(D); for (int i = 0; i < D; i++) { @@ -280,10 +275,10 @@ std::vector equifactor(int const N, int const D, int const randseed = 0) { } return products; -} // equipartition +} // equifactor -#ifdef ENABLE_DEBUG +#ifndef NDEBUG void print_sets(std::vector> sets) { int nsets = sets.size(); int minprod = std::numeric_limits::max(); diff --git a/wonton/support/lsfits.h b/wonton/support/lsfits.h index e5a27e0..338ecc9 100644 --- a/wonton/support/lsfits.h +++ b/wonton/support/lsfits.h @@ -26,6 +26,100 @@ Please see the license file at the root of this repository, or at: namespace Wonton { +/** + * @brief Build the matrices used in the least square + * equation (A^T.A).X = (A^T.F) from the stencil. + * + * Since the stencil is the same for multiple field variables, + * it should be called only once. It will also be used to compute + * the right hand side of the equation. Here the first point + * of the stencil is assumed to be the point where the gradient + * must be computed. + * + * @tparam dim: the spatial dimension of the problem. + * @param stencil: the stencil point coordinates. + * @param invert: invert (A^T.A) or not. + * @return (A^T.A) or its inverse, and A^T. + */ +template +std::vector build_gradient_stencil_matrices(std::vector> const& stencil, + bool invert = false) { + + // the first point is the reference point where we are trying + // to compute the gradient; so the number of rows is size - 1. + int const num = stencil.size() - 1; + + // each row of A contains the components of the vector from + // coords[0] to the candidate point coords[i] being used + // in the least squares approximation (x_i - x_0). + Matrix A(num, dim); + for (int i = 0; i < num; ++i) { + for (int j = 0; j < dim; ++j) { + A[i][j] = stencil[i+1][j] - stencil[0][j]; + } + } + + Matrix AT = A.transpose(); + return { invert ? (AT * A).inverse() : (AT * A), AT }; +} + +/** + * @brief Build the right hand side (A^T.F) of the least square + * equation (A^T.A).X = (A^T.F) from field values. + * + * It should be called for each field variable using the cached + * transposed matrix A^T used in the least square equation. + * + * @tparam dim: the spatial dimension of the problem. + * @param AT: the cached transposed matrix A^T. + * @param values: the field value at each stencil point. + * @return the vector (A^T.F) + */ +template +Vector build_right_hand_side(Matrix const& AT, std::vector const& values) { + + // the first value is that of the reference point where we are trying + // to compute the gradient; so the number of components is size - 1. + int const num = values.size() - 1; + + // Each entry/row of F contains the difference between the + // function value at the candidate point and the function value + // at the point where we are computing (f-f_0) + std::vector F(num); + for (int i = 0; i < num; ++i) { + F[i] = values[i+1] - values[0]; + } + + return Vector(AT * F); +} + +/** + * @brief Compute the gradient using the cached inverse of the + * (A^T.A) matrix and the given right hand side (A^T.F). + * + * @tparam dim: the spatial dimension of the problem. + * @tparam CoordSys: the current coordinate system. + * @param ATA_inv: the cached inverse of the least square matrix. + * @param AT: the cached transposed matrix A^T. + * @param values: the field value at each stencil point. + * @return the gradient at the current point. + */ +template +Vector ls_gradient(Matrix const& ATA_inv, + Matrix const& AT, + std::vector const& values) { + + CoordSys::template verify_coordinate_system(); + + // compute the right hand side (A^T.F) from stencil values + Vector ATF = build_right_hand_side(AT, values); + + // compute gradient + Vector gradient = ATA_inv * ATF; + + return gradient; +} + /*! @brief Compute least squares gradient from set of values @param[in] coords Vector of coordinates at which values are given @@ -36,6 +130,9 @@ namespace Wonton { and the first value is assumed to the value at this reference point This operator does not know anything about a mesh. + It solves the algebraic equation using an optimized LAPACK kernel + if available, or the inverse method if not. It is meant to be used + only when remapping a single field variable. */ @@ -45,62 +142,26 @@ Vector ls_gradient(std::vector> const & coords, CoordSys::template verify_coordinate_system(); - Point coord0 = coords[0]; - - double val0 = vals[0]; - - // There are nvals but the first is the reference point where we - // are trying to compute the gradient; so the matrix sizes etc - // will only be nvals-1 - - int nvals = vals.size(); - - // Each row of A contains the components of the vector from - // coord0 to the candidate point being used in the Least Squares - // approximation (X_i-X_0). - - Matrix A(nvals-1, D); - for (int i = 0; i < nvals-1; ++i) { - for (int j = 0; j < D; ++j) - A[i][j] = coords[i+1][j]-coord0[j]; - } - - - // A is a matrix of size nvals-1 by D (where D is the space - // dimension). So transpose(A)*A is D by D - - Matrix AT = A.transpose(); - - Matrix ATA = AT*A; - - // Each entry/row of F contains the difference between the - // function value at the candidate point and the function value - // at the point where we are computing (f-f_0) - - std::vector F(nvals-1); - for (int i = 0; i < nvals-1; ++i) - F[i] = vals[i+1]-val0; - - // F is a vector of nvals. So transpose(A)*F is vector of D - // (where D is the space dimension) + // construct the least square equation components + auto M = build_gradient_stencil_matrices(coords, false); + Vector ATF = build_right_hand_side(M[1], vals); - Vector ATF = Vector(AT*F); + // solve it + Matrix& ATA = M[0]; - // Inverse of ATA +#ifdef WONTON_HAS_LAPACKE + Matrix B(D, 1); + for (int i = 0; i < D; ++i) { B[i][0] = ATF[i]; } - Matrix ATAinv = ATA.inverse(); - - // Gradient of length D - - auto gradient = ATAinv*ATF; - - // Corrections for curvilinear coordinates - - CoordSys::modify_gradient(gradient, coord0); - - // Return + // use lapack solver for symmetric positive definite matrices + Vector gradient( ATA.solve(B, "lapack-posv").data() ); +#else + // use the QR decomposition method + Vector gradient; + solve(ATA, ATF, gradient); +#endif - return ATAinv*ATF; + return gradient; } diff --git a/wonton/support/test/test_polytope_2.cc b/wonton/support/test/test_polytope_2.cc index 3e067dc..abbda77 100644 --- a/wonton/support/test/test_polytope_2.cc +++ b/wonton/support/test/test_polytope_2.cc @@ -72,11 +72,11 @@ TEST(Polytope, Mesh2D) { Wonton::Point<2>(0.0, 1.0) }; - std::vector ncv_poly_moments = { 0.75, 0.375, 21.0/72.0 }; + std::vector ncv_poly_moments = { 0.75, 0.375, 21.0/72.0, 25.0/96.0, 14.0/96.0, 15.0/96.0}; Wonton::Polytope<2> ncv_poly(ncv_poly_points); - // Verify moments - auto poly_moments = ncv_poly.moments(); - for (int im = 0; im < 3; im++) + // Verify moments upto order 2 + auto poly_moments = ncv_poly.moments(2); + for (int im = 0; im < 6; im++) ASSERT_NEAR(ncv_poly_moments[im], poly_moments[im], 1.0e-15); }