From 8741434c02ac47165522f04f61beed95a6a625d6 Mon Sep 17 00:00:00 2001 From: acrovato <a.crovato@uliege.be> Date: Tue, 1 Feb 2022 13:44:36 +0100 Subject: [PATCH] Rm waves, add amfe. Fix format. --- .gitlab-ci.yml | 100 +++++++++++++++++----------------------- .gitmodules | 3 ++ CMakeLists.txt | 10 ++++ README.md | 4 +- ext/CMakeLists.txt | 1 + ext/amfe | 1 + fpm/_src/CMakeLists.txt | 24 +++++----- fpm/src/CMakeLists.txt | 8 +--- fpm/src/fBody.cpp | 8 +++- fpm/src/fBuilder.cpp | 67 +++++++++++++-------------- fpm/src/fBuilder.h | 3 +- fpm/src/fField.cpp | 8 ++++ fpm/src/fField.h | 2 +- fpm/src/fProblem.cpp | 87 +++++++++++++++++----------------- fpm/src/fSolver.cpp | 39 +++++++--------- fpm/src/fWake.cpp | 6 ++- fpm/tests/n12t.py | 8 ++-- run.py | 36 +++++++++++++++ run_fpm.py | 26 ----------- 19 files changed, 225 insertions(+), 216 deletions(-) create mode 100644 .gitmodules create mode 100644 ext/CMakeLists.txt create mode 160000 ext/amfe create mode 100755 run.py delete mode 100755 run_fpm.py diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index eaec999..cb1c837 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -11,80 +11,66 @@ default: - mn2l variables: - GIT_CLONE_PATH: /builds/$CI_PROJECT_PATH/FPM_BUILD/fpm # cannot work in /builds/$CI_PROJECT_PATH + GIT_SUBMODULE_STRATEGY: recursive + GIT_STRATEGY: clone # workaround full clone for each pipeline (https://gitlab.com/gitlab-org/gitlab-runner/-/issues/26993) + GIT_LFS_SKIP_SMUDGE: 1 # do not pull LFS stages: - build -# - fmt_dox + - test build: <<: *global_tag_def stage: build script: - - printenv | sort - - cd .. - # get and build waves - - rm -rf waves - - wget -q https://gitlab.uliege.be/am-dept/waves/-/archive/master/waves-master.tar.bz2 - - tar xf waves-master.tar.bz2 - - rm waves-master.tar.bz2 - - mv waves-master waves - - cd waves + - git submodule init + - git submodule update + - rm -rf build workspace - mkdir build - cd build - - cmake -Wno-dev -C ../CMake/disable-trilinos.cmake .. + - cmake -Wno-dev .. - make -j $(nproc) - # check code format + artifacts: + paths: + - build/ + expire_in: 1 day + +format: + <<: *global_tag_def + stage: build + script: - clang-format --version # we use clang-format-10 exclusively - - cd ../../fpm - - ${CI_PROJECT_DIR}/../waves/scripts/format_code.py + - ./ext/amfe/scripts/format_code.py - mkdir -p patches - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi - # build fpm - - mkdir build - - cd build - - cmake -Wno-dev .. # -DCMAKE_PREFIX_PATH=${CI_PROJECT_DIR}/../waves (handled by default) - - make -j $(nproc) - # build documentation - - make dox - # test - - ctest -j $(nproc) --output-on-failure #--verbose - - mv ${CI_PROJECT_DIR}/../waves/scripts/format_code.py . # ulgy way to keep a script we need later... - #timeout: 10 hours # will be available in 12.3 artifacts: paths: - - build/ - patches/ expire_in: 1 day + when: on_failure + dependencies: + - build + allow_failure: true -#format: -# <<: *global_tag_def -# stage: fmt_dox -# script: -# - clang-format --version # we use clang-format-10 exclusively -# - ./build/format_code.py -# - mkdir -p patches -# - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi -# artifacts: -# paths: -# - patches/ -# expire_in: 1 day -# when: on_failure -# dependencies: -# - build_test -# allow_failure: true -# -# -#doxygen: -# <<: *global_tag_def -# stage: fmt_dox -# script: -# - cd build -# - make dox -# artifacts: -# paths: -# - build/doxygen/ -# expire_in: 1 day -# dependencies: -# - build_test +doxygen: + <<: *global_tag_def + stage: test + script: + - cd build + - make dox + artifacts: + paths: + - build/doxygen/ + expire_in: 1 week + dependencies: + - build +ctest: + <<: *global_tag_def + stage: test + script: + - cd build + - ctest --output-on-failure -j $(nproc) + #timeout: 10 hours # will be available in 12.3 + dependencies: + - build diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..e84a099 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "ext/amfe"] + path = ext/amfe + url = ../amfe.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 9df63d9..0258585 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -116,9 +116,19 @@ ELSE() MESSAGE("Doxygen need to be installed to generate the doxygen documentation") ENDIF() +# -- Define (for SWIG to detect definitions) +INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}) # to find "amfe_def.h" + # -- CTest ENABLE_TESTING() # -- Modules +ADD_SUBDIRECTORY( ext ) ADD_SUBDIRECTORY( fpm ) +# -- Final +MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}") +MESSAGE(STATUS "* SYSTEM NAME=\"${CMAKE_SYSTEM_NAME}\"") +MESSAGE(STATUS "* CXX COMPILER: ${CMAKE_CXX_COMPILER_ID}") +MESSAGE(STATUS "* CXX STANDARD: ${CMAKE_CXX_STANDARD}") +MESSAGE(STATUS "* BUILD TYPE: ${CMAKE_BUILD_TYPE}") diff --git a/README.md b/README.md index 2212ff5..1650e53 100644 --- a/README.md +++ b/README.md @@ -10,12 +10,12 @@ Set of python/C++ modules solving the linear potential equation and correcting i - [x] [eigen3](http://eigen.tuxfamily.org/) support for linear algebra ## Build -fpm is linked against [waves](https://gitlab.uliege.be/am-dept/waves). Once all the required packages and waves are installed: +fpm depends on [amfe](https://gitlab.uliege.be/am-dept/amfe). Once all the required packages are installed: ```bash git clone git@gitlab.uliege.be:am-dept/fpm.git # get the code mkdir build && cd build -cmake -DCMAKE_PREFIX_PATH=/path/to/waves .. # configure and reference waves package +cmake .. # configure make # compile ctest # test ``` diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt new file mode 100644 index 0000000..fddb82f --- /dev/null +++ b/ext/CMakeLists.txt @@ -0,0 +1 @@ +ADD_SUBDIRECTORY( amfe) \ No newline at end of file diff --git a/ext/amfe b/ext/amfe new file mode 160000 index 0000000..054ba00 --- /dev/null +++ b/ext/amfe @@ -0,0 +1 @@ +Subproject commit 054ba00c327e4f2745ff9c559832247f90c3cec9 diff --git a/fpm/_src/CMakeLists.txt b/fpm/_src/CMakeLists.txt index a137601..bd2af2e 100644 --- a/fpm/_src/CMakeLists.txt +++ b/fpm/_src/CMakeLists.txt @@ -22,16 +22,12 @@ FILE(GLOB ISRCS *.i) SET(CMAKE_SWIG_FLAGS "") SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON) -FOREACH(dir ${WAVES_INCLUDE_DIRS}) - LIST(APPEND WAVES_SINCLUDE_DIRS "-I${dir}") -ENDFOREACH() -FOREACH(dir ${WAVES_SWIG_DIRS}) - LIST(APPEND WAVES_SINCLUDE_DIRS "-I${dir}") -ENDFOREACH() - -SET(SWINCFLAGS +SET(SWINCFLAGS -I${PROJECT_SOURCE_DIR}/fpm/src -${WAVES_SINCLUDE_DIRS} +-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/src +-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src +-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/src +-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src ) SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}") @@ -43,11 +39,13 @@ endif() MACRO_DebugPostfix(_fpmw) TARGET_INCLUDE_DIRECTORIES(_fpmw PRIVATE ${PROJECT_SOURCE_DIR}/fpm/_src - ${WAVES_SWIG_DIRS} - ${PYTHON_INCLUDE_PATH} + ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src + ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src + ${PYTHON_INCLUDE_PATH} ) -SWIG_LINK_LIBRARIES(fpmw - fpm ${WAVES_LIBRARIES} ${PYTHON_LIBRARIES} +SWIG_LINK_LIBRARIES(fpmw + fpm tbox fwk ${PYTHON_LIBRARIES} ) + diff --git a/fpm/src/CMakeLists.txt b/fpm/src/CMakeLists.txt index 15f283e..6664a4f 100644 --- a/fpm/src/CMakeLists.txt +++ b/fpm/src/CMakeLists.txt @@ -26,12 +26,6 @@ FIND_PACKAGE(EIGEN 3.3.4 REQUIRED) TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${EIGEN_INCLUDE_DIRS}) TARGET_COMPILE_DEFINITIONS(fpm PUBLIC EIGEN_DONT_PARALLELIZE) -# -- WAVES (if no path is provided, assume that it is located next to fpm) -IF(NOT DEFINED CMAKE_PREFIX_PATH) - SET(CMAKE_PREFIX_PATH "${PROJECT_SOURCE_DIR}/../waves/build") -ENDIF() -FIND_PACKAGE(WAVES REQUIRED) -TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${WAVES_INCLUDE_DIRS}) -TARGET_LINK_LIBRARIES(fpm ${WAVES_LIBRARIES}) +TARGET_LINK_LIBRARIES(fpm tbox) SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)") diff --git a/fpm/src/fBody.cpp b/fpm/src/fBody.cpp index 2b9bf3e..398e3c0 100644 --- a/fpm/src/fBody.cpp +++ b/fpm/src/fBody.cpp @@ -30,6 +30,10 @@ using namespace fpm; Body::Body(std::shared_ptr<MshData> _msh, std::string const &id, std::vector<std::string> const &teIds, double xF) : Group(_msh, id), Cl(0), Cd(0), Cs(0), Cm(0) { + // Initialize element memory + for (auto e : tag->elems) + e->initValues(false); + // If wakes are requested, check if they are already available, otherwise create them bool hasChanged = false; size_t j = 0; @@ -70,7 +74,7 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id, err << "fpm::Body: zero or negative length wake requested for " << *tag << "!\n"; throw std::runtime_error(err.str()); } - Node *nN = new Node(msh->nodes.back()->no + 1, Eigen::Vector3d(xF, n->pos(1), n->pos(2))); + Node *nN = new Node(msh->nodes.back()->no + 1, Eigen::Vector3d(xF, n->pos(1), n->pos(2)), msh->nodes.back()->row + 1); wkNodes.push_back(nN); msh->nodes.push_back(nN); } @@ -85,7 +89,7 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id, std::map<Node *, Node *> teMap; for (auto n : teNodes) { - Node *nN = new Node(msh->nodes.back()->no + 1, n->pos); + Node *nN = new Node(msh->nodes.back()->no + 1, n->pos, msh->nodes.back()->row + 1); teMap[n] = nN; msh->nodes.push_back(nN); } diff --git a/fpm/src/fBuilder.cpp b/fpm/src/fBuilder.cpp index fe5a573..e349f55 100644 --- a/fpm/src/fBuilder.cpp +++ b/fpm/src/fBuilder.cpp @@ -23,7 +23,6 @@ #include "wTag.h" #include "wElement.h" #include "wNode.h" -#include "wMem.h" #define PI 3.14159 @@ -80,7 +79,7 @@ void Builder::run() Eigen::MatrixXd trsfCoordSrc(3, nVertices); // Computations related to source panels - glob2loc = Glob2loc(false, ej); // Transformation matrix + glob2loc = Glob2loc(false, ej); // Transformation matrix trsfCoordSrc = TrsfCoordSrc(false, ej, glob2loc); // Local coordinates of the source panel // body panels @@ -94,8 +93,8 @@ void Builder::run() Eigen::Vector2d mutau; // Computations related to target panels - trsfCoordTgt = TrsfCoordTgt(false, ei, ej, glob2loc); // Local coordinates of the target panel - distance = Distance(false, ei, ej, glob2loc); // Distances between source and target panels + trsfCoordTgt = TrsfCoordTgt(false, ei, ej, glob2loc); // Local coordinates of the target panel + distance = Distance(false, ei, ej, glob2loc); // Distances between source and target panels mutau = MuTau(false, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC // Assembly of the matrices @@ -119,7 +118,7 @@ void Builder::run() Eigen::MatrixXd trsfCoordSrc(3, nVertices); // Computations related to source panels - glob2loc = Glob2loc(true, ej); // Transformation matrix + glob2loc = Glob2loc(true, ej); // Transformation matrix trsfCoordSrc = TrsfCoordSrc(true, ej, glob2loc); // Local coordinates of the source panel // body panels @@ -133,8 +132,8 @@ void Builder::run() Eigen::Vector2d mutau; // Computations related to target panels - trsfCoordTgt = TrsfCoordTgt(true, ei, ej, glob2loc); // Local coordinates of the target panel - distance = Distance(true, ei, ej, glob2loc); // Distances between source and target panels + trsfCoordTgt = TrsfCoordTgt(true, ei, ej, glob2loc); // Local coordinates of the target panel + distance = Distance(true, ei, ej, glob2loc); // Distances between source and target panels mutau = MuTau(true, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC // Assembly of the matrices @@ -159,7 +158,7 @@ void Builder::run() Eigen::MatrixXd trsfCoordSrc(3, nVertices); // Computations related to source panels - glob2loc = Glob2loc(false, ew); // Transformation matrix + glob2loc = Glob2loc(false, ew); // Transformation matrix trsfCoordSrc = TrsfCoordSrc(false, ew, glob2loc); // Local coordinates of the source panel // wake panels @@ -173,8 +172,8 @@ void Builder::run() Eigen::Vector2d wmutau; // Computations related to target panels - trsfCoordTgt = TrsfCoordTgt(false, ei, ew, glob2loc); // Local coordinates of the target panel - distance = Distance(false, ei, ew, glob2loc); // Distances between source and target panels + trsfCoordTgt = TrsfCoordTgt(false, ei, ew, glob2loc); // Local coordinates of the target panel + distance = Distance(false, ei, ew, glob2loc); // Distances between source and target panels wmutau = MuTau(false, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC // Assembly of the matrices @@ -201,7 +200,7 @@ void Builder::run() Eigen::MatrixXd trsfCoordSrc(3, nVertices); // Computations related to source panels - glob2loc = Glob2loc(true, ew); // Transformation matrix + glob2loc = Glob2loc(true, ew); // Transformation matrix trsfCoordSrc = TrsfCoordSrc(true, ew, glob2loc); // Local coordinates of the source panel // wake panels @@ -215,8 +214,8 @@ void Builder::run() Eigen::Vector2d wmutau; // Computations related to target panels - trsfCoordTgt = TrsfCoordTgt(true, ei, ew, glob2loc); // Local coordinates of the target panel - distance = Distance(true, ei, ew, glob2loc); // Distances between source and target panels + trsfCoordTgt = TrsfCoordTgt(true, ei, ew, glob2loc); // Local coordinates of the target panel + distance = Distance(true, ei, ew, glob2loc); // Distances between source and target panels wmutau = MuTau(true, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC // Assembly of the matrices @@ -249,11 +248,11 @@ void Builder::run() */ Eigen::Matrix3d Builder::Glob2loc(bool symy, Element *ej) { - Eigen::Matrix3d glob2loc; // Transformation matrix - Eigen::Vector3d n = ej->normal(); // Unit normal vector of the source panel + Eigen::Matrix3d glob2loc; // Transformation matrix + Eigen::Vector3d n = ej->getNormal(); // Unit normal vector of the source panel // Computation of the unit longitudinal vector - Source panel - Eigen::Vector3d cg = ej->getVMem().getCg(); + Eigen::Vector3d cg = ej->getCg(); Eigen::Vector3d x0 = ej->nodes[0]->pos; // Coordinates of node 1 Eigen::Vector3d l = (x0 - cg).normalized(); @@ -265,7 +264,7 @@ Eigen::Matrix3d Builder::Glob2loc(bool symy, Element *ej) glob2loc.row(1) = p; glob2loc.row(2) = n; - if(symy) + if (symy) { glob2loc.row(1) = -glob2loc.row(1); glob2loc.col(1) = -glob2loc.col(1); @@ -287,7 +286,7 @@ Eigen::MatrixXd Builder::TrsfCoordSrc(bool symy, Element *ej, Eigen::Matrix3d co for (int j = 0; j < nVertices; ++j) coords(i, j) = ej->nodes[(nVertices - 1) - j]->pos[i]; // Global coordinates of the nodes of the source panel - if(symy) + if (symy) coords.row(1) = -coords.row(1); for (int i = 0; i < nVertices; ++i) @@ -302,8 +301,8 @@ Eigen::MatrixXd Builder::TrsfCoordSrc(bool symy, Element *ej, Eigen::Matrix3d co Eigen::Vector3d Builder::TrsfCoordTgt(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc) { Eigen::Vector3d trsfCoordTgt; - Eigen::Vector3d cg = ei->getVMem().getCg(); // Global coordinates of the CG of the target panel - trsfCoordTgt = glob2loc * cg; // Local coordinates of the target panel + Eigen::Vector3d cg = ei->getCg(); // Global coordinates of the CG of the target panel + trsfCoordTgt = glob2loc * cg; // Local coordinates of the target panel return trsfCoordTgt; } @@ -314,13 +313,13 @@ Eigen::Vector3d Builder::TrsfCoordTgt(bool symy, Element *ei, Element *ej, Eigen Eigen::Vector3d Builder::Distance(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc) { Eigen::Vector3d distance; // Distance between the target and source panels - Eigen::Vector3d cgSrc = ej->getVMem().getCg(); - Eigen::Vector3d cgTgt = ei->getVMem().getCg(); + Eigen::Vector3d cgSrc = ej->getCg(); + Eigen::Vector3d cgTgt = ei->getCg(); for (int i = 0; i < 3; ++i) distance(i) = cgTgt(i) - cgSrc(i); - if(symy) + if (symy) distance(1) = cgTgt(1) + 2 * cgSrc(1); distance = glob2loc * distance; @@ -347,24 +346,24 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri // Coefficients for (int i = 1; i <= nVertices; ++i) { - e(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + distance(2) * distance(2); - h(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)); - r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2)); + e(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + distance(2) * distance(2); + h(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)); + r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) + distance(2) * distance(2)); } for (int i = 1; i <= nVertices; ++i) { int j = i % nVertices + 1; - F(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(i - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(i - 1); - G(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(j - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(j - 1); + F(i - 1) = (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * e(i - 1) - (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * h(i - 1); + G(i - 1) = (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * e(j - 1) - (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * h(j - 1); } // Doublet term - Influence of a panel on another panel for (int i = 1; i <= nVertices; ++i) { int j = i % nVertices + 1; - mu = mu + atan2(distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)), - distance(2) * distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1)); + mu = mu + atan2(distance(2) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)), + distance(2) * distance(2) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1)); } mu = -1 / (4 * PI) * mu; } @@ -375,8 +374,8 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri { int j = i % nVertices + 1; - d(i - 1) = sqrt((trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) + (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1))); - r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2)); + d(i - 1) = sqrt((trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) + (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1))); + r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) + distance(2) * distance(2)); } if (ei->no == ej->no && !symy) @@ -386,7 +385,7 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri { int j = i % nVertices + 1; - tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1))); + tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1))); } tau = -1 / (4 * PI) * tau; } @@ -397,7 +396,7 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri { int j = i % nVertices + 1; - tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1))); + tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1))); } tau = -1 / (4 * PI) * tau - distance(2) * mu; diff --git a/fpm/src/fBuilder.h b/fpm/src/fBuilder.h index 0e721a8..10069df 100644 --- a/fpm/src/fBuilder.h +++ b/fpm/src/fBuilder.h @@ -36,8 +36,7 @@ namespace fpm */ class FPM_API Builder : public fwk::wSharedObject { - -std::map<tbox::Element *, int> rows; ///< map linking an element to a matrix cell + std::map<tbox::Element *, int> rows; ///< map linking an element to a matrix cell public: std::shared_ptr<Problem> pbl; ///< problem definition diff --git a/fpm/src/fField.cpp b/fpm/src/fField.cpp index 37f76c1..b0f06a0 100644 --- a/fpm/src/fField.cpp +++ b/fpm/src/fField.cpp @@ -15,11 +15,19 @@ */ #include "fField.h" +#include "wElement.h" #include "wTag.h" using namespace tbox; using namespace fpm; +Field::Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name) : Group(_msh, name) +{ + // Initialize element memory + for (auto e : tag->elems) + e->initValues(true); +} + void Field::write(std::ostream &out) const { out << "fpm::Field " << *tag << std::endl; diff --git a/fpm/src/fField.h b/fpm/src/fField.h index 6ef5bb8..360c422 100644 --- a/fpm/src/fField.h +++ b/fpm/src/fField.h @@ -30,7 +30,7 @@ namespace fpm class FPM_API Field : public tbox::Group { public: - Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name) : Group(_msh, name) {} + Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name); ~Field() { std::cout << "~Field()\n"; } #ifndef SWIG diff --git a/fpm/src/fProblem.cpp b/fpm/src/fProblem.cpp index 5d1b840..f240905 100644 --- a/fpm/src/fProblem.cpp +++ b/fpm/src/fProblem.cpp @@ -22,7 +22,6 @@ #include "wTag.h" #include "wElement.h" #include "wNode.h" -#include "wMem.h" #include "wCache.h" using namespace tbox; using namespace fpm; @@ -57,69 +56,69 @@ void Problem::add(std::shared_ptr<Body> b) Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, double const &tau) { Eigen::Vector3d U(0, 0, 0); - size_t nVertices = e.nodes.size(); // Number of nodes per panel - size_t nGP = e.getVCache().getVGauss().getN(); // Number of Gauss points on the element - Eigen::MatrixXd coords(3,nVertices); + size_t nVertices = e.nodes.size(); // Number of nodes per panel + size_t nGP = e.getVCache().getVGauss().getN(); // Number of Gauss points on the element + Eigen::MatrixXd coords(3, nVertices); for (int i = 0; i < 3; ++i) for (int j = 0; j < nVertices; ++j) - coords(i,j) = e.nodes[(nVertices - 1) - j]->pos[i]; // Coordinates of the panel - + coords(i, j) = e.nodes[(nVertices - 1) - j]->pos[i]; // Coordinates of the panel + // Import of the shape functions at the Gauss points - Eigen::MatrixXd ff(nVertices,nGP); - for(int i = 0; i < nVertices; ++i) - for(int j = 0; j < nGP; ++j) - ff(i,j) = e.getVCache().getSf(j)[i]; + Eigen::MatrixXd ff(nVertices, nGP); + for (int i = 0; i < nVertices; ++i) + for (int j = 0; j < nGP; ++j) + ff(i, j) = e.getVCache().getSf(j)[i]; // Import of the derivatives of the shape functions at the Gauss points - Eigen::MatrixXd dff(2 * nGP,nVertices); - for(int i = 0; i < nGP; ++i) - for(int j = 0; j < 2; ++j) + Eigen::MatrixXd dff(2 * nGP, nVertices); + for (int i = 0; i < nGP; ++i) + for (int j = 0; j < 2; ++j) dff.row(j + (2 * i)) = e.getVCache().getDsf(i).row(j); - + // Computation of the Jacobian matrix - xi-derivative (first row of the Jacobian) - Eigen::MatrixXd dXi(3,nGP); - for(int i = 0; i < nGP; ++i) + Eigen::MatrixXd dXi(3, nGP); + for (int i = 0; i < nGP; ++i) { dXi.col(i) << 0, 0, 0; - for(int j = 0; j < nVertices; ++j) + for (int j = 0; j < nVertices; ++j) { - dXi(0,i) += dff((2 * i),j) * coords(0,j); - dXi(1,i) += dff((2 * i),j) * coords(1,j); - dXi(2,i) += dff((2 * i),j) * coords(2,j); + dXi(0, i) += dff((2 * i), j) * coords(0, j); + dXi(1, i) += dff((2 * i), j) * coords(1, j); + dXi(2, i) += dff((2 * i), j) * coords(2, j); } } // Computation of the Jacobian matrix - eta-derivative (second row of the Jacobian) - Eigen::MatrixXd dEta(3,nGP); - for(int i = 0; i < nGP; ++i) + Eigen::MatrixXd dEta(3, nGP); + for (int i = 0; i < nGP; ++i) { dEta.col(i) << 0, 0, 0; - for(int j = 0; j < nVertices; ++j) + for (int j = 0; j < nVertices; ++j) { - dEta(0,i) += dff(1 + (2 * i),j) * coords(0,j); - dEta(1,i) += dff(1 + (2 * i),j) * coords(1,j); - dEta(2,i) += dff(1 + (2 * i),j) * coords(2,j); + dEta(0, i) += dff(1 + (2 * i), j) * coords(0, j); + dEta(1, i) += dff(1 + (2 * i), j) * coords(1, j); + dEta(2, i) += dff(1 + (2 * i), j) * coords(2, j); } } // Computation of the Jacobian matrix - cross-product (third row of the Jacobian) - Eigen::MatrixXd dZeta(3,nGP); + Eigen::MatrixXd dZeta(3, nGP); - for(int i = 0; i < nGP; ++i) + for (int i = 0; i < nGP; ++i) { Eigen::Vector3d dxi = dXi.col(i); Eigen::Vector3d deta = dEta.col(i); Eigen::Vector3d dzeta; dzeta = (dxi.cross(deta)).normalized(); - dZeta.col(i) = dzeta; + dZeta.col(i) = dzeta; } // Assembly of the Jacobian and partial computation of the velocity // Done for each Gauss point -> For n Gauss points, there are n values of the velocity -> Average - Eigen::MatrixXd v(3,nGP); - for(int i = 0; i < nGP; ++i) + Eigen::MatrixXd v(3, nGP); + for (int i = 0; i < nGP; ++i) { // Creation of the Jacobian Eigen::Matrix3d J0; @@ -130,9 +129,9 @@ Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, doub // Inversion of the Jacobian Eigen::Matrix3d Jinv = J0.inverse(); - // Removal of the third column - Eigen::MatrixXd J(3,2); - for(int j = 0; j < 2; ++j) + // Removal of the third column + Eigen::MatrixXd J(3, 2); + for (int j = 0; j < 2; ++j) { J.col(j) = Jinv.col(j); } @@ -145,28 +144,28 @@ Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, doub } // Computation of the velocity - Eigen::MatrixXd dff0(2,4); + Eigen::MatrixXd dff0(2, 4); dff0.row(0) = dff.row(2 * i); dff0.row(1) = dff.row(1 + (2 * i)); v.col(i) = J * dff0 * mu_element; } // Average of the velocity - for(int i = 0; i < nGP; ++i) + for (int i = 0; i < nGP; ++i) { U += v.col(i); } - U = U/nGP; + U = U / nGP; // Import of the unit normal vector - Target panel - Eigen::Vector3d n = e.normal(); + Eigen::Vector3d n = e.getNormal(); // Contribution of the normal velocity - U += (-tau)*n; + U += (-tau) * n; // Contribution of the freestream velocity U += Ui(); - + return U; } @@ -221,17 +220,17 @@ void Problem::updateMem() // Update volume Jacobian if (field) for (auto e : field->tag->elems) - e->getVMem().update(true); + e->update(); // Update surface Jacobian for (auto s : bodies) { for (auto e : s->tag->elems) { - e->getVMem().update(false); + e->update(); } - for(auto w : s-> wakes) + for (auto w : s->wakes) for (auto e : w->tag->elems) - e->getVMem().update(false); + e->update(); } } diff --git a/fpm/src/fSolver.cpp b/fpm/src/fSolver.cpp index 6e49acc..36328b2 100644 --- a/fpm/src/fSolver.cpp +++ b/fpm/src/fSolver.cpp @@ -22,7 +22,6 @@ #include "wMshData.h" #include "wNode.h" #include "wElement.h" -#include "wMem.h" #include "wTag.h" #include "wResults.h" #include "wMshExport.h" @@ -36,20 +35,14 @@ using namespace fpm; /** * @brief Initialize the solver and perform sanity checks - * @authors Adrien Crovato & Axel Dechamps */ Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0) { // Say hi std::cout << std::endl; std::cout << "*******************************************************************************" << std::endl; - std::cout << "** \\_/ **" << std::endl; - std::cout << "** \\_______O(_)O_______/ **" << std::endl; - std::cout << "** **" << std::endl; - std::cout << "*******************************************************************************" << std::endl; - std::cout << "** Hi! My name is FPM v0.1-2007 **" << std::endl; - std::cout << "** Adrien Crovato & Romain Boman **" << std::endl; - std::cout << "** ULiege 2020-2021 **" << std::endl; + std::cout << "** Hi! My name is FPM v1.0-2201 **" << std::endl; + std::cout << "** ULiege 2020-2022 **" << std::endl; std::cout << "*******************************************************************************" << std::endl << std::endl; @@ -89,8 +82,8 @@ void Solver::run() { for (int i = 0; i < body->tag->elems.size(); ++i, ++ig) { - tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal()); - taue[body->tag->elems[i]->no-1] = aic->pbl->Ui().dot(body->tag->elems[i]->normal()); + tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->getNormal()); + taue[body->tag->elems[i]->no - 1] = aic->pbl->Ui().dot(body->tag->elems[i]->getNormal()); } } @@ -148,7 +141,7 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, Eigen::VectorXd phie = -aic->A * mu - aic->B * tau - aic->C * sigma; mue.resize(mu.size()); - for(size_t i = 0 ; i < mue.size() ; ++i) + for (size_t i = 0; i < mue.size(); ++i) mue[i] = mu(i); auto pbl = aic->pbl; @@ -173,10 +166,10 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, // average at nodes for (auto n : body->nodes) { - for (auto e : body->neMap.at(n)) + for (auto e : body->neMap.at(n)) { - vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); - muNodes[n->no - 1] += mapMu.at(e) / body->neMap.at(n).size(); + vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); + muNodes[n->no - 1] += mapMu.at(e) / body->neMap.at(n).size(); } } } @@ -190,10 +183,10 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, int ig = 0; for (auto e : body->tag->elems) { - U[e->no - 1] = pbl->U(*e, muNodes, tau(ig)); // velocity - rho[e->no - 1] = pbl->Rho(U[e->no - 1]); // density - M[e->no - 1] = pbl->Mach(U[e->no - 1]); // Mach number - Cp[e->no - 1] = pbl->Cp(U[e->no - 1]); // pressure coefficient + U[e->no - 1] = pbl->U(*e, muNodes, tau(ig)); // velocity + rho[e->no - 1] = pbl->Rho(U[e->no - 1]); // density + M[e->no - 1] = pbl->Mach(U[e->no - 1]); // Mach number + Cp[e->no - 1] = pbl->Cp(U[e->no - 1]); // pressure coefficient ig += 1; } } @@ -219,7 +212,7 @@ void Solver::computeLoad() body->cLoadZ[i] = 0; for (auto e : body->neMap.at(body->nodes[i])) { - Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVMem().getVol() * e->normal() / e->nodes.size(); + Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVol() * e->getNormal() / e->nodes.size(); body->cLoadX[i] += cLoad(0); body->cLoadY[i] += cLoad(1); body->cLoadZ[i] += cLoad(2); @@ -229,11 +222,11 @@ void Solver::computeLoad() // Compute aerodynamic load coefficients (i.e. Load / pressure / surface) // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D) body->Cm = 0; - Eigen::Vector3d Cf(0, 0, 0); + Eigen::Vector3d Cf(0, 0, 0); for (auto e : body->tag->elems) { - Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR; // lever arm - Eigen::Vector3d cf = Cp[e->no - 1] * e->getVMem().getVol() * e->normal(); // force coefficient + Eigen::Vector3d l = e->getCg() - pbl->xR; // lever arm + Eigen::Vector3d cf = Cp[e->no - 1] * e->getVol() * e->getNormal(); // force coefficient Cf -= cf; body->Cm -= (l(2) * cf(0) - l(0) * cf(2)) / (pbl->sR * pbl->cR); // moment is positive along y (3D) and negative along z (2D) => positive nose-up } diff --git a/fpm/src/fWake.cpp b/fpm/src/fWake.cpp index 0b40761..f8ba042 100644 --- a/fpm/src/fWake.cpp +++ b/fpm/src/fWake.cpp @@ -25,6 +25,10 @@ using namespace fpm; Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &wTag) : Group(_msh, name) { + // Initialize element memory + for (auto e : tag->elems) + e->initValues(false); + // Create set of TE edges std::unordered_set<Edge *, EquEdge, EquEdge> teEdges; for (auto e : tag->elems) @@ -50,7 +54,7 @@ Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &w if (it != teEdges.end()) { // check normal direction - if (e->normal().dot((*it)->el0->normal()) > 0) + if (e->getNormal().dot((*it)->el0->getNormal()) > 0) (*it)->el1 = e; else (*it)->el2 = e; diff --git a/fpm/tests/n12t.py b/fpm/tests/n12t.py index fa1d843..64a73a8 100644 --- a/fpm/tests/n12t.py +++ b/fpm/tests/n12t.py @@ -27,7 +27,7 @@ def main(): p['sym'] = 0 p['spn'] = 5 # true span if sym=0, half span if sym=1 p['tspn'] = 3 # true span if sym=0, half span if sym=1 - p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'tSpn' : p['tspn'], 'tNc' : 100, 'tNs' : 10} + p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'tSpn' : p['tspn'], 'tNc' : 100, 'tNs' : 5} # flow parameters p['aoa'] = 3 * math.pi / 180 p['aos'] = 0 @@ -39,10 +39,10 @@ def main(): print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) tests = CTests() if(p['aoa'] == 3 * math.pi / 180 and p['mach'] == 0): - tests.add(CTest('CL', solver.Cl, 0.273, 5e-2)) - tests.add(CTest('CD', solver.Cd, 0.004, 5e-2)) + tests.add(CTest('CL', solver.Cl, 0.241, 5e-2)) + tests.add(CTest('CD', solver.Cd, 0.0045, 5e-2)) tests.add(CTest('CS', solver.Cs, -0.000, 5e-2)) - tests.add(CTest('CM', solver.Cm, -0.277, 5e-2)) + tests.add(CTest('CM', solver.Cm, -0.113, 5e-2)) tests.run() # eof diff --git a/run.py b/run.py new file mode 100755 index 0000000..86fb30f --- /dev/null +++ b/run.py @@ -0,0 +1,36 @@ +#!/usr/bin/env python3 +# -*- coding: utf8 -*- +# test encoding: à -é-è-ô-ï-€ + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Calls fwk.wutils.run to execute a script as if dartflo was installed + +def main(): + import os.path, sys + # adds fwk/tbox to the python path + thisdir = os.path.split(os.path.abspath(__file__))[0] + fwkdir = os.path.abspath(os.path.join(thisdir, 'ext', 'amfe')) + if not os.path.isdir(fwkdir): + raise Exception('fpm/ext/amfe not found!\n') + sys.path.append(fwkdir) + # adds "." to the pythonpath + sys.path.append(thisdir) + + import fwk.wutils as wu + wu.run() + +if __name__ == "__main__": + main() diff --git a/run_fpm.py b/run_fpm.py deleted file mode 100755 index a63a0ec..0000000 --- a/run_fpm.py +++ /dev/null @@ -1,26 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf8 -*- -# test encoding: à -é-è-ô-ï-€ -# -# runs a test as if it was installed -# - fixes the python path in a dev environment -# - creates a workspace folder - -if __name__ == "__main__": - import sys - import os - - thisdir = os.path.split(os.path.abspath(__file__))[0] - - # look for fwk - fwkdir = os.path.abspath(os.path.join(thisdir, '..', 'waves')) - if not os.path.isdir(fwkdir): - raise Exception( - "'waves' not found - clone it next to fpm (from https://gitlab.uliege.be/am-dept/waves)") - sys.path.append(fwkdir) - - # adds "." to the pythonpath (after waves so that waves is found first) - sys.path.append(thisdir) - - import run # maybe we should put the code of run in fwk.... - run.main(thisdir) -- GitLab