From 035acd9b27999c7d18a268bb21bd38fcf6be02d3 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Thu, 12 Jan 2023 17:27:20 +0100
Subject: [PATCH] Add steady sdpm and tests. Simplify and Refactor mesh
management class. Update docker image and clang. Update credits. Format.
---
.gitlab-ci.yml | 4 +-
CREDITS.md | 1 +
scripts/format_code.py | 3 +-
sdpm/_src/sdpmw.i | 64 ++++---
sdpm/models/oneraM6.geo | 327 +++++++++++++++++++++++++++++++++++
sdpm/src/sdpm.h | 10 +-
sdpm/src/sdpmBody.cpp | 150 ++++++++++++++++
sdpm/src/sdpmBody.h | 65 +++++++
sdpm/src/sdpmBuilder.cpp | 139 +++++++++++++++
sdpm/src/sdpmBuilder.h | 40 +++++
sdpm/src/sdpmEdge.h | 91 ++++++++++
sdpm/src/sdpmElement.cpp | 28 +--
sdpm/src/sdpmElement.h | 22 +--
sdpm/src/sdpmElement.inl | 56 ------
sdpm/src/sdpmGmshExport.cpp | 119 +++++++------
sdpm/src/sdpmGmshExport.h | 7 +-
sdpm/src/sdpmGmshImport.cpp | 22 +--
sdpm/src/sdpmGmshImport.h | 4 +-
sdpm/src/sdpmGroup.cpp | 40 +++++
sdpm/src/sdpmGroup.h | 48 ++++++
sdpm/src/sdpmLine2.cpp | 36 +---
sdpm/src/sdpmLine2.h | 11 +-
sdpm/src/sdpmMesh.cpp | 22 +--
sdpm/src/sdpmMesh.h | 26 ++-
sdpm/src/sdpmNode.cpp | 3 +-
sdpm/src/sdpmNode.h | 2 +-
sdpm/src/sdpmProblem.cpp | 112 ++++++++++++
sdpm/src/sdpmProblem.h | 86 +++++++++
sdpm/src/sdpmQuad4.cpp | 69 ++++----
sdpm/src/sdpmQuad4.h | 13 +-
sdpm/src/sdpmSolver.cpp | 336 ++++++++++++++++++++++++++++++++++++
sdpm/src/sdpmSolver.h | 76 ++++++++
sdpm/src/sdpmTimers.h | 4 +-
sdpm/src/sdpmUnitTest.cpp | 23 +--
sdpm/src/sdpmUnitTest.h | 2 +-
sdpm/src/sdpmWake.cpp | 156 +++++++++++++++++
sdpm/src/sdpmWake.h | 60 +++++++
sdpm/src/sdpm_extras.h | 2 +-
sdpm/tests/gmshio.py | 22 +--
sdpm/tests/m6.py | 68 ++++++++
sdpm/tests/n12.py | 70 ++++++++
sdpm/tests/n12tail.py | 77 +++++++++
42 files changed, 2194 insertions(+), 322 deletions(-)
create mode 100644 sdpm/models/oneraM6.geo
create mode 100644 sdpm/src/sdpmBody.cpp
create mode 100644 sdpm/src/sdpmBody.h
create mode 100644 sdpm/src/sdpmBuilder.cpp
create mode 100644 sdpm/src/sdpmBuilder.h
create mode 100644 sdpm/src/sdpmEdge.h
delete mode 100644 sdpm/src/sdpmElement.inl
create mode 100644 sdpm/src/sdpmGroup.cpp
create mode 100644 sdpm/src/sdpmGroup.h
create mode 100644 sdpm/src/sdpmProblem.cpp
create mode 100644 sdpm/src/sdpmProblem.h
create mode 100644 sdpm/src/sdpmSolver.cpp
create mode 100644 sdpm/src/sdpmSolver.h
create mode 100644 sdpm/src/sdpmWake.cpp
create mode 100644 sdpm/src/sdpmWake.h
create mode 100644 sdpm/tests/m6.py
create mode 100644 sdpm/tests/n12.py
create mode 100644 sdpm/tests/n12tail.py
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8d447ab..570f357 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,5 +1,5 @@
default:
- image: rboman/waves-py3:2020.3
+ image: rboman/waves-py3:2022.0
before_script:
- source /opt/intel/mkl/bin/mklvars.sh intel64
- source /opt/intel/tbb/bin/tbbvars.sh intel64
@@ -22,7 +22,7 @@ format:
<<: *global_tag_def
stage: build
script:
- - clang-format --version # we use clang-format-10 exclusively
+ - clang-format --version
- ./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
diff --git a/CREDITS.md b/CREDITS.md
index 77ab805..8e5959f 100644
--- a/CREDITS.md
+++ b/CREDITS.md
@@ -6,3 +6,4 @@ Direct and indirect contributions have been made to SDPM by the following people
- Grigorios Dimitriadis
- Romain Boman
- Luc Papeleux
+- Axel Dechamps
diff --git a/scripts/format_code.py b/scripts/format_code.py
index 908f80e..36b03f9 100755
--- a/scripts/format_code.py
+++ b/scripts/format_code.py
@@ -40,13 +40,12 @@ def main():
encs = {}
for f in all_files(os.getcwd(), patterns='*.cpp;*.c;*.h;*.hpp'):
# print(f)
- cmd = ['clang-format-10', "-style=file", "-i", f]
+ cmd = ['clang-format', "-style=file", "-i", f]
retcode = subprocess.call(cmd)
if retcode != 0:
print(f'ERROR: retcode = {retcode}')
break
-
if __name__ == "__main__":
# print('running format_code.py...')
main()
diff --git a/sdpm/_src/sdpmw.i b/sdpm/_src/sdpmw.i
index 7cdaac3..8a384a2 100644
--- a/sdpm/_src/sdpmw.i
+++ b/sdpm/_src/sdpmw.i
@@ -14,7 +14,7 @@
* limitations under the License.
*/
-// fwk.i: SWIG input file for python interface
+// SWIG input file for python interface
%feature("autodoc","1");
@@ -33,10 +33,15 @@ threads="1"
#include "sdpmLine2.h"
#include "sdpmQuad4.h"
#include "sdpmTag.h"
+#include "sdpmGroup.h"
#include "sdpmGmshImport.h"
#include "sdpmGmshExport.h"
#include "sdpmUnitTest.h"
+#include "sdpmBody.h"
+#include "sdpmProblem.h"
+#include "sdpmSolver.h"
+
#include "sdpmCppBuf2Py.h"
#include <sstream>
@@ -51,34 +56,34 @@ threads="1"
// from: http://swig.10945.n7.nabble.com/Trapping-Swig-DirectorException-td6013.html
// le code suivant permet de voir la call stack dans les appels C++ => python
%include "exception.i"
-%{
- static void handle_exception(void) {
- try {
- throw;
- } catch (std::exception &e) {
- std::stringstream txt;
+%{
+ static void handle_exception(void) {
+ try {
+ throw;
+ } catch (std::exception &e) {
+ std::stringstream txt;
txt << e.what(); // << ' ' << typeid(e).name();
- PyErr_SetString(PyExc_RuntimeError, e.what());
- }
- catch(...)
+ PyErr_SetString(PyExc_RuntimeError, e.what());
+ }
+ catch(...)
{
PyErr_SetString(PyExc_RuntimeError, "Unknown C++ Runtime Error");
- }
- }
-%}
-%exception {
- try {
- $action
- } catch (...) {
- // Note that if a director method failed, the Python error indicator
- // already contains full details of the exception, and it will be
- // reraised when we go to SWIG_fail; so no need to convert the C++
- // exception back to a Python one
- if (!PyErr_Occurred()) {
- handle_exception();
- }
- SWIG_fail;
- }
+ }
+ }
+%}
+%exception {
+ try {
+ $action
+ } catch (...) {
+ // Note that if a director method failed, the Python error indicator
+ // already contains full details of the exception, and it will be
+ // reraised when we go to SWIG_fail; so no need to convert the C++
+ // exception back to a Python one
+ if (!PyErr_Occurred()) {
+ handle_exception();
+ }
+ SWIG_fail;
+ }
}
// SDPM
@@ -92,11 +97,16 @@ threads="1"
%include "sdpmLine2.h"
%include "sdpmQuad4.h"
%include "sdpmTag.h"
+%include "sdpmGroup.h"
%include "sdpmGmshImport.h"
%include "sdpmGmshExport.h"
%include "sdpmUnitTest.h"
-%warnfilter(401); //Nothing known about base class 'std::basic_streambuf< char >'
+%include "sdpmBody.h"
+%include "sdpmProblem.h"
+%include "sdpmSolver.h"
+
+%warnfilter(401); // Nothing known about base class 'std::basic_streambuf< char >'
%include "sdpmCppBuf2Py.h"
%warnfilter(+401);
diff --git a/sdpm/models/oneraM6.geo b/sdpm/models/oneraM6.geo
new file mode 100644
index 0000000..418d940
--- /dev/null
+++ b/sdpm/models/oneraM6.geo
@@ -0,0 +1,327 @@
+/* Onera M6 wing */
+
+// Mesh parameters
+// domain and mesh
+DefineConstant[ wNc = { 50, Name "Number of panel along chord" }
+ wNs = { 10, Name "Number of panel along span" }
+ bC = { 0.05, Name "Progression along chord" }
+ pS = { 1.0, Name "Progression along span" } ];
+
+//// GEOMETRY
+/// Points
+// Airfoil 1: oneraM6, 143 points
+Point(1) = {0.805900,0.000000,0.000000};
+Point(2) = {0.804129,0.000000,0.000232};
+Point(3) = {0.802038,0.000000,0.000505};
+Point(4) = {0.799569,0.000000,0.000828};
+Point(5) = {0.796652,0.000000,0.001210};
+Point(6) = {0.793208,0.000000,0.001661};
+Point(7) = {0.789139,0.000000,0.002193};
+Point(8) = {0.784331,0.000000,0.002822};
+Point(9) = {0.778649,0.000000,0.003565};
+Point(10) = {0.771932,0.000000,0.004444};
+Point(11) = {0.763989,0.000000,0.005483};
+Point(12) = {0.754592,0.000000,0.006710};
+Point(13) = {0.743470,0.000000,0.008151};
+Point(14) = {0.730299,0.000000,0.009828};
+Point(15) = {0.714691,0.000000,0.011690};
+Point(16) = {0.696182,0.000000,0.013774};
+Point(17) = {0.677546,0.000000,0.015783};
+Point(18) = {0.658809,0.000000,0.017717};
+Point(19) = {0.639970,0.000000,0.019586};
+Point(20) = {0.621026,0.000000,0.021397};
+Point(21) = {0.601980,0.000000,0.023159};
+Point(22) = {0.582826,0.000000,0.024875};
+Point(23) = {0.563568,0.000000,0.026547};
+Point(24) = {0.544202,0.000000,0.028170};
+Point(25) = {0.524729,0.000000,0.029737};
+Point(26) = {0.505147,0.000000,0.031238};
+Point(27) = {0.485455,0.000000,0.032658};
+Point(28) = {0.465652,0.000000,0.033984};
+Point(29) = {0.445738,0.000000,0.035197};
+Point(30) = {0.425711,0.000000,0.036282};
+Point(31) = {0.405570,0.000000,0.037225};
+Point(32) = {0.385317,0.000000,0.038011};
+Point(33) = {0.364947,0.000000,0.038631};
+Point(34) = {0.344460,0.000000,0.039073};
+Point(35) = {0.323856,0.000000,0.039344};
+Point(36) = {0.303135,0.000000,0.039432};
+Point(37) = {0.282293,0.000000,0.039343};
+Point(38) = {0.261331,0.000000,0.039078};
+Point(39) = {0.240248,0.000000,0.038642};
+Point(40) = {0.219042,0.000000,0.038037};
+Point(41) = {0.197712,0.000000,0.037261};
+Point(42) = {0.176258,0.000000,0.036306};
+Point(43) = {0.154679,0.000000,0.035154};
+Point(44) = {0.132972,0.000000,0.033774};
+Point(45) = {0.111131,0.000000,0.032117};
+Point(46) = {0.092662,0.000000,0.030457};
+Point(47) = {0.077073,0.000000,0.028830};
+Point(48) = {0.063937,0.000000,0.027269};
+Point(49) = {0.052891,0.000000,0.025800};
+Point(50) = {0.043619,0.000000,0.024441};
+Point(51) = {0.035851,0.000000,0.023203};
+Point(52) = {0.029356,0.000000,0.022027};
+Point(53) = {0.023936,0.000000,0.020812};
+Point(54) = {0.019428,0.000000,0.019503};
+Point(55) = {0.015684,0.000000,0.018096};
+Point(56) = {0.012586,0.000000,0.016619};
+Point(57) = {0.010032,0.000000,0.015114};
+Point(58) = {0.007931,0.000000,0.013618};
+Point(59) = {0.006214,0.000000,0.012165};
+Point(60) = {0.004815,0.000000,0.010776};
+Point(61) = {0.003683,0.000000,0.009463};
+Point(62) = {0.002775,0.000000,0.008233};
+Point(63) = {0.002050,0.000000,0.007089};
+Point(64) = {0.001480,0.000000,0.006027};
+Point(65) = {0.001037,0.000000,0.005045};
+Point(66) = {0.000698,0.000000,0.004138};
+Point(67) = {0.000444,0.000000,0.003301};
+Point(68) = {0.000260,0.000000,0.002529};
+Point(69) = {0.000135,0.000000,0.001818};
+Point(70) = {0.000056,0.000000,0.001162};
+Point(71) = {0.000013,0.000000,0.000557};
+Point(72) = {0.000000,0.000000,0.000000};
+Point(73) = {0.000013,0.000000,-0.000557};
+Point(74) = {0.000056,0.000000,-0.001162};
+Point(75) = {0.000135,0.000000,-0.001818};
+Point(76) = {0.000260,0.000000,-0.002529};
+Point(77) = {0.000444,0.000000,-0.003301};
+Point(78) = {0.000698,0.000000,-0.004138};
+Point(79) = {0.001037,0.000000,-0.005045};
+Point(80) = {0.001480,0.000000,-0.006027};
+Point(81) = {0.002050,0.000000,-0.007089};
+Point(82) = {0.002775,0.000000,-0.008233};
+Point(83) = {0.003683,0.000000,-0.009463};
+Point(84) = {0.004815,0.000000,-0.010776};
+Point(85) = {0.006214,0.000000,-0.012165};
+Point(86) = {0.007931,0.000000,-0.013618};
+Point(87) = {0.010032,0.000000,-0.015114};
+Point(88) = {0.012586,0.000000,-0.016619};
+Point(89) = {0.015684,0.000000,-0.018096};
+Point(90) = {0.019428,0.000000,-0.019503};
+Point(91) = {0.023936,0.000000,-0.020812};
+Point(92) = {0.029356,0.000000,-0.022027};
+Point(93) = {0.035851,0.000000,-0.023203};
+Point(94) = {0.043619,0.000000,-0.024441};
+Point(95) = {0.052891,0.000000,-0.025800};
+Point(96) = {0.063937,0.000000,-0.027269};
+Point(97) = {0.077073,0.000000,-0.028830};
+Point(98) = {0.092662,0.000000,-0.030457};
+Point(99) = {0.111131,0.000000,-0.032117};
+Point(100) = {0.132972,0.000000,-0.033774};
+Point(101) = {0.154679,0.000000,-0.035154};
+Point(102) = {0.176258,0.000000,-0.036306};
+Point(103) = {0.197712,0.000000,-0.037261};
+Point(104) = {0.219042,0.000000,-0.038037};
+Point(105) = {0.240248,0.000000,-0.038642};
+Point(106) = {0.261331,0.000000,-0.039078};
+Point(107) = {0.282293,0.000000,-0.039343};
+Point(108) = {0.303135,0.000000,-0.039432};
+Point(109) = {0.323856,0.000000,-0.039344};
+Point(110) = {0.344460,0.000000,-0.039073};
+Point(111) = {0.364947,0.000000,-0.038631};
+Point(112) = {0.385317,0.000000,-0.038011};
+Point(113) = {0.405570,0.000000,-0.037225};
+Point(114) = {0.425711,0.000000,-0.036282};
+Point(115) = {0.445738,0.000000,-0.035197};
+Point(116) = {0.465652,0.000000,-0.033984};
+Point(117) = {0.485455,0.000000,-0.032658};
+Point(118) = {0.505147,0.000000,-0.031238};
+Point(119) = {0.524729,0.000000,-0.029737};
+Point(120) = {0.544202,0.000000,-0.028170};
+Point(121) = {0.563568,0.000000,-0.026547};
+Point(122) = {0.582826,0.000000,-0.024875};
+Point(123) = {0.601980,0.000000,-0.023159};
+Point(124) = {0.621026,0.000000,-0.021397};
+Point(125) = {0.639970,0.000000,-0.019586};
+Point(126) = {0.658809,0.000000,-0.017717};
+Point(127) = {0.677546,0.000000,-0.015783};
+Point(128) = {0.696182,0.000000,-0.013774};
+Point(129) = {0.714691,0.000000,-0.011690};
+Point(130) = {0.730299,0.000000,-0.009828};
+Point(131) = {0.743470,0.000000,-0.008151};
+Point(132) = {0.754592,0.000000,-0.006710};
+Point(133) = {0.763989,0.000000,-0.005483};
+Point(134) = {0.771932,0.000000,-0.004444};
+Point(135) = {0.778649,0.000000,-0.003565};
+Point(136) = {0.784331,0.000000,-0.002822};
+Point(137) = {0.789139,0.000000,-0.002193};
+Point(138) = {0.793208,0.000000,-0.001661};
+Point(139) = {0.796652,0.000000,-0.001210};
+Point(140) = {0.799569,0.000000,-0.000828};
+Point(141) = {0.802038,0.000000,-0.000505};
+Point(142) = {0.804129,0.000000,-0.000232};
+// Airfoil 2: oneraM6, 143 points
+Point(144) = {1.143427,1.196000,0.000000};
+Point(145) = {1.142432,1.196000,0.000130};
+Point(146) = {1.141256,1.196000,0.000284};
+Point(147) = {1.139869,1.196000,0.000466};
+Point(148) = {1.138230,1.196000,0.000680};
+Point(149) = {1.136294,1.196000,0.000933};
+Point(150) = {1.134007,1.196000,0.001232};
+Point(151) = {1.131305,1.196000,0.001586};
+Point(152) = {1.128112,1.196000,0.002004};
+Point(153) = {1.124337,1.196000,0.002498};
+Point(154) = {1.119873,1.196000,0.003082};
+Point(155) = {1.114592,1.196000,0.003771};
+Point(156) = {1.108341,1.196000,0.004581};
+Point(157) = {1.100939,1.196000,0.005523};
+Point(158) = {1.092167,1.196000,0.006570};
+Point(159) = {1.081765,1.196000,0.007741};
+Point(160) = {1.071292,1.196000,0.008870};
+Point(161) = {1.060762,1.196000,0.009957};
+Point(162) = {1.050174,1.196000,0.011007};
+Point(163) = {1.039528,1.196000,0.012025};
+Point(164) = {1.028824,1.196000,0.013015};
+Point(165) = {1.018059,1.196000,0.013980};
+Point(166) = {1.007236,1.196000,0.014919};
+Point(167) = {0.996353,1.196000,0.015831};
+Point(168) = {0.985409,1.196000,0.016712};
+Point(169) = {0.974403,1.196000,0.017556};
+Point(170) = {0.963336,1.196000,0.018354};
+Point(171) = {0.952208,1.196000,0.019099};
+Point(172) = {0.941016,1.196000,0.019781};
+Point(173) = {0.929760,1.196000,0.020391};
+Point(174) = {0.918441,1.196000,0.020920};
+Point(175) = {0.907059,1.196000,0.021362};
+Point(176) = {0.895611,1.196000,0.021711};
+Point(177) = {0.884097,1.196000,0.021959};
+Point(178) = {0.872518,1.196000,0.022111};
+Point(179) = {0.860873,1.196000,0.022161};
+Point(180) = {0.849160,1.196000,0.022111};
+Point(181) = {0.837379,1.196000,0.021962};
+Point(182) = {0.825530,1.196000,0.021717};
+Point(183) = {0.813612,1.196000,0.021377};
+Point(184) = {0.801625,1.196000,0.020941};
+Point(185) = {0.789568,1.196000,0.020404};
+Point(186) = {0.777440,1.196000,0.019757};
+Point(187) = {0.765241,1.196000,0.018981};
+Point(188) = {0.752966,1.196000,0.018050};
+Point(189) = {0.742587,1.196000,0.017117};
+Point(190) = {0.733826,1.196000,0.016203};
+Point(191) = {0.726444,1.196000,0.015325};
+Point(192) = {0.720236,1.196000,0.014500};
+Point(193) = {0.715025,1.196000,0.013736};
+Point(194) = {0.710659,1.196000,0.013040};
+Point(195) = {0.707009,1.196000,0.012379};
+Point(196) = {0.703963,1.196000,0.011696};
+Point(197) = {0.701429,1.196000,0.010961};
+Point(198) = {0.699325,1.196000,0.010170};
+Point(199) = {0.697584,1.196000,0.009340};
+Point(200) = {0.696149,1.196000,0.008494};
+Point(201) = {0.694968,1.196000,0.007654};
+Point(202) = {0.694003,1.196000,0.006837};
+Point(203) = {0.693217,1.196000,0.006056};
+Point(204) = {0.692581,1.196000,0.005318};
+Point(205) = {0.692070,1.196000,0.004627};
+Point(206) = {0.691663,1.196000,0.003984};
+Point(207) = {0.691343,1.196000,0.003387};
+Point(208) = {0.691094,1.196000,0.002835};
+Point(209) = {0.690903,1.196000,0.002325};
+Point(210) = {0.690760,1.196000,0.001855};
+Point(211) = {0.690657,1.196000,0.001421};
+Point(212) = {0.690587,1.196000,0.001022};
+Point(213) = {0.690542,1.196000,0.000653};
+Point(214) = {0.690518,1.196000,0.000313};
+Point(215) = {0.690511,1.196000,0.000000};
+Point(216) = {0.690518,1.196000,-0.000313};
+Point(217) = {0.690542,1.196000,-0.000653};
+Point(218) = {0.690587,1.196000,-0.001022};
+Point(219) = {0.690657,1.196000,-0.001421};
+Point(220) = {0.690760,1.196000,-0.001855};
+Point(221) = {0.690903,1.196000,-0.002325};
+Point(222) = {0.691094,1.196000,-0.002835};
+Point(223) = {0.691343,1.196000,-0.003387};
+Point(224) = {0.691663,1.196000,-0.003984};
+Point(225) = {0.692070,1.196000,-0.004627};
+Point(226) = {0.692581,1.196000,-0.005318};
+Point(227) = {0.693217,1.196000,-0.006056};
+Point(228) = {0.694003,1.196000,-0.006837};
+Point(229) = {0.694968,1.196000,-0.007654};
+Point(230) = {0.696149,1.196000,-0.008494};
+Point(231) = {0.697584,1.196000,-0.009340};
+Point(232) = {0.699325,1.196000,-0.010170};
+Point(233) = {0.701429,1.196000,-0.010961};
+Point(234) = {0.703963,1.196000,-0.011696};
+Point(235) = {0.707009,1.196000,-0.012379};
+Point(236) = {0.710659,1.196000,-0.013040};
+Point(237) = {0.715025,1.196000,-0.013736};
+Point(238) = {0.720236,1.196000,-0.014500};
+Point(239) = {0.726444,1.196000,-0.015325};
+Point(240) = {0.733826,1.196000,-0.016203};
+Point(241) = {0.742587,1.196000,-0.017117};
+Point(242) = {0.752966,1.196000,-0.018050};
+Point(243) = {0.765241,1.196000,-0.018981};
+Point(244) = {0.777440,1.196000,-0.019757};
+Point(245) = {0.789568,1.196000,-0.020404};
+Point(246) = {0.801625,1.196000,-0.020941};
+Point(247) = {0.813612,1.196000,-0.021377};
+Point(248) = {0.825530,1.196000,-0.021717};
+Point(249) = {0.837379,1.196000,-0.021962};
+Point(250) = {0.849160,1.196000,-0.022111};
+Point(251) = {0.860873,1.196000,-0.022161};
+Point(252) = {0.872518,1.196000,-0.022111};
+Point(253) = {0.884097,1.196000,-0.021959};
+Point(254) = {0.895611,1.196000,-0.021711};
+Point(255) = {0.907059,1.196000,-0.021362};
+Point(256) = {0.918441,1.196000,-0.020920};
+Point(257) = {0.929760,1.196000,-0.020391};
+Point(258) = {0.941016,1.196000,-0.019781};
+Point(259) = {0.952208,1.196000,-0.019099};
+Point(260) = {0.963336,1.196000,-0.018354};
+Point(261) = {0.974403,1.196000,-0.017556};
+Point(262) = {0.985409,1.196000,-0.016712};
+Point(263) = {0.996353,1.196000,-0.015831};
+Point(264) = {1.007236,1.196000,-0.014919};
+Point(265) = {1.018059,1.196000,-0.013980};
+Point(266) = {1.028824,1.196000,-0.013015};
+Point(267) = {1.039528,1.196000,-0.012025};
+Point(268) = {1.050174,1.196000,-0.011007};
+Point(269) = {1.060762,1.196000,-0.009957};
+Point(270) = {1.071292,1.196000,-0.008870};
+Point(271) = {1.081765,1.196000,-0.007741};
+Point(272) = {1.092167,1.196000,-0.006570};
+Point(273) = {1.100939,1.196000,-0.005523};
+Point(274) = {1.108341,1.196000,-0.004581};
+Point(275) = {1.114592,1.196000,-0.003771};
+Point(276) = {1.119873,1.196000,-0.003082};
+Point(277) = {1.124337,1.196000,-0.002498};
+Point(278) = {1.128112,1.196000,-0.002004};
+Point(279) = {1.131305,1.196000,-0.001586};
+Point(280) = {1.134007,1.196000,-0.001232};
+Point(281) = {1.136294,1.196000,-0.000933};
+Point(282) = {1.138230,1.196000,-0.000680};
+Point(283) = {1.139869,1.196000,-0.000466};
+Point(284) = {1.141256,1.196000,-0.000284};
+Point(285) = {1.142432,1.196000,-0.000130};
+
+/// Lines
+// Airfoil 1:
+Spline(1) = {1:72};
+Spline(2) = {72:142,1};
+// Airfoil 2:
+Spline(3) = {144:215};
+Spline(4) = {215:285,144};
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,144};
+Line(62) = {72,215};
+
+/// Line loops & Surfaces
+// Wing 1:
+Line Loop(11) = {1,62,-3,-61};
+Line Loop(12) = {2,61,-4,-62};
+Surface(11) = {-11};
+Surface(12) = {-12};
+
+//// MESH
+Transfinite Line{1,2,3,4} = wNc + 1 Using Bump bC;
+Transfinite Line{61,62} = wNs + 1 Using Progression pS;
+Transfinite Surface{11,12};
+Recombine Surface{11,12};
+
+//// PHYSICAL GROUPS
+// Trailing edge
+Physical Line("wTe") = {61};
+// Wing:
+Physical Surface("wing") = {11,12};
diff --git a/sdpm/src/sdpm.h b/sdpm/src/sdpm.h
index 40c8b53..4fe47a5 100644
--- a/sdpm/src/sdpm.h
+++ b/sdpm/src/sdpm.h
@@ -36,6 +36,7 @@
using sdpmDouble = double;
using sdpmMatrixXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, Eigen::Dynamic>;
using sdpmVectorXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, 1>;
+using sdpmMatrix3d = Eigen::Matrix<sdpmDouble, 3, 3>;
using sdpmVector3d = Eigen::Matrix<sdpmDouble, 3, 1>;
namespace sdpm
@@ -51,10 +52,17 @@ class Element;
class Line2;
class Quad4;
class Tag;
+class Group;
class Result;
class GmshImport;
class GmshExport;
// Panel method
+class Problem;
+class Body;
+class Wake;
+class Edge;
+class Builder;
+class Solver;
} // namespace sdpm
-#endif //SDPM_H
+#endif // SDPM_H
diff --git a/sdpm/src/sdpmBody.cpp b/sdpm/src/sdpmBody.cpp
new file mode 100644
index 0000000..6d85032
--- /dev/null
+++ b/sdpm/src/sdpmBody.cpp
@@ -0,0 +1,150 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpm.h"
+#include "sdpmBody.h"
+#include "sdpmMesh.h"
+#include "sdpmElement.h"
+#include "sdpmQuad4.h"
+#include "sdpmTag.h"
+#include "sdpmGmshExport.h"
+#include <iostream>
+using namespace sdpm;
+
+Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double xF) : Group(msh, name), cl(0), cd(0), cs(0), cm(0)
+{
+ // Create wake
+ std::cout << "Creating wake... " << std::flush;
+ std::string wkName = _tag->getName() + "Wake";
+ // get TE nodes and sort them according to their y-coordinate (ascending)
+ Group te(_msh, teName);
+ std::vector<Node *> teNodes;
+ for (auto e : te.getElements())
+ for (auto n : e->getNodes())
+ teNodes.push_back(n);
+ std::sort(teNodes.begin(), teNodes.end(), [](Node *a, Node *b) -> bool
+ { return a->getCoords()(1) < b->getCoords()(1); });
+ teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
+ // create wake geometry if it does not already exist
+ std::map<std::string, Tag *> tags = _msh.getTags();
+ auto it = tags.find(wkName);
+ if (it == tags.end())
+ {
+ // create tag and add it to the mesh
+ Tag *tagp = new Tag(tags.rbegin()->second->getId() + 1, wkName, 2);
+ _msh.addTag(wkName, tagp);
+ // translate TE nodes (along x-coordinate)
+ std::vector<Node *> wkNodes;
+ for (auto n : teNodes)
+ {
+ sdpmVector3d pos = n->getCoords();
+ // sanity check
+ if (xF <= pos(0))
+ {
+ std::stringstream err;
+ err << "sdpm::Body: zero or negative length wake requested for " << *_tag << "!\n";
+ throw std::runtime_error(err.str());
+ }
+ Node *nN = new Node(_msh.getNodes().back()->getId() + 1, sdpmVector3d(xF, pos(1), pos(2)));
+ _msh.addNode(nN);
+ wkNodes.push_back(nN);
+ }
+ // create wake elements
+ for (size_t i = 0; i < wkNodes.size() - 1; ++i)
+ {
+ std::vector<Node *> qnodes = {teNodes[i], wkNodes[i], wkNodes[i + 1], teNodes[i + 1]};
+ _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
+ }
+ // duplicate TE nodes
+ std::map<Node *, Node *> teMap;
+ for (auto n : teNodes)
+ {
+ Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
+ _msh.addNode(nN);
+ teMap[n] = nN;
+ }
+ // create wake
+ _wake = new Wake(_msh, wkName, _tag);
+ // swap nodes
+ for (auto p : _wake->getMap())
+ {
+ Element *lwEl = p.second.second;
+ std::vector<Node *> lwNods = lwEl->getNodes();
+ for (size_t i = 0; i < lwNods.size(); ++i)
+ {
+ try
+ {
+ lwEl->setNode(i, teMap.at(lwNods[i]));
+ }
+ catch (const std::out_of_range &)
+ {
+ // std::cout << lwe->nodes[i]->no << "not found in map!\n";
+ }
+ }
+ }
+ // get body nodes (now that they are duplicated and swaped at lower TE)
+ getUniqueNodes();
+ // print and save new mesh
+ std::cout << *_wake << "created. " << std::flush;
+ GmshExport writer(_msh);
+ writer.save();
+ }
+ else
+ {
+ // get body nodes
+ getUniqueNodes();
+ // map the lower TE nodes to their upper nodes
+ std::map<Node *, Node *> iTeMap;
+ for (auto ten : teNodes)
+ for (auto wn : _nodes)
+ if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
+ iTeMap[wn] = ten;
+ // create wake and print
+ _wake = new Wake(_msh, wkName, _tag, iTeMap);
+ std::cout << *_wake << "already existing, nothing done." << std::endl;
+ }
+
+ // Associate each node to its adjacent elements
+ for (auto e : _tag->getElements())
+ for (auto n : e->getNodes())
+ _neMap[n].push_back(e);
+
+ // Size load vectors
+ nLoads.resize(_nodes.size());
+}
+
+Body::~Body()
+{
+ delete _wake;
+ std::cout << "~sdpm::Body()\n";
+}
+
+/**
+ * @brief Return a vector of unique nodes belonging to the tag
+ */
+void Body::getUniqueNodes()
+{
+ for (auto e : _tag->getElements())
+ for (auto n : e->getNodes())
+ _nodes.push_back(n);
+ std::sort(_nodes.begin(), _nodes.end());
+ _nodes.erase(std::unique(_nodes.begin(), _nodes.end()), _nodes.end());
+}
+
+void Body::write(std::ostream &out) const
+{
+ out << "sdpm::Body " << *_tag << "with " << *_wake;
+}
diff --git a/sdpm/src/sdpmBody.h b/sdpm/src/sdpmBody.h
new file mode 100644
index 0000000..c3fa2c8
--- /dev/null
+++ b/sdpm/src/sdpmBody.h
@@ -0,0 +1,65 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMBODY_H
+#define SDPMBODY_H
+
+#include "sdpm.h"
+#include "sdpmGroup.h"
+#include "sdpmWake.h"
+#include <map>
+
+namespace sdpm
+{
+
+/**
+ * @brief Manage a lifting body immersed in the fluid
+ * @authors Adrien Crovato
+ */
+class SDPM_API Body : public sdpm::Group
+{
+ std::vector<Node *> _nodes; ///< nodes of the surface
+ std::map<Node *, std::vector<Element *>> _neMap; ///< map between nodes and adjacent elements
+ Wake *_wake; ///< wake attached to the lifting body
+
+public:
+ sdpmDouble cl; ///< lift coefficient
+ sdpmDouble cd; ///< drag coefficient
+ sdpmDouble cs; ///< sideforce coefficient
+ sdpmDouble cm; ///< pitch moment coefficient (positive nose-up)
+ std::vector<sdpmVector3d> nLoads; ///< nodal aerodynamic load (normalized by dynamic pressure)
+
+ Body(Mesh &msh, std::string const &name, std::string const &teName, double xF);
+ ~Body();
+
+private:
+ void getUniqueNodes();
+
+public:
+ std::vector<Node *> const &getNodes() const { return _nodes; }
+#ifndef SWIG
+ std::map<Node *, std::vector<Element *>> const &getMap() const
+ {
+ return _neMap;
+ }
+ Wake const &getWake() const { return *_wake; }
+ virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace sdpm
+
+#endif // SDPMBODY_H
diff --git a/sdpm/src/sdpmBuilder.cpp b/sdpm/src/sdpmBuilder.cpp
new file mode 100644
index 0000000..709718f
--- /dev/null
+++ b/sdpm/src/sdpmBuilder.cpp
@@ -0,0 +1,139 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpmBuilder.h"
+#include "sdpmElement.h"
+#include "sdpmNode.h"
+
+using namespace sdpm;
+
+/**
+ * @brief Compute local frame of source panel
+ */
+sdpmMatrix3d Builder::computeTransfoMat(bool symy, Element const *ej)
+{
+ // Compute tangents and normal
+ sdpmVector3d u = (ej->getNodes()[0]->getCoords() - ej->getCg()).normalized();
+ sdpmVector3d v_ = ej->getNodes()[1]->getCoords() - ej->getCg();
+ sdpmVector3d v = (v_ - (u.dot(v_)) * u).normalized();
+ sdpmVector3d n = ej->getNormal();
+
+ // Assemble into transformation matrix
+ sdpmMatrix3d g2l;
+ g2l.row(0) = u;
+ g2l.row(1) = v;
+ g2l.row(2) = n;
+ if (symy)
+ {
+ g2l.row(1) = -g2l.row(1);
+ g2l.col(1) = -g2l.col(1);
+ }
+ return g2l;
+}
+
+/**
+ * @brief Gather the coordinates of source panel
+ */
+sdpmMatrixXd Builder::gatherCoords(bool symy, Element const *ej)
+{
+ // Get nodes
+ std::vector<Node *> nods = ej->getNodes();
+ size_t nV = nods.size();
+
+ // Get coordinates
+ Eigen::MatrixXd coords(3, nV);
+ for (int j = 0; j < nV; ++j)
+ coords.col(j) = nods[nV - 1 - j]->getCoords();
+ if (symy)
+ coords.row(1) = -coords.row(1);
+ return coords;
+}
+
+/**
+ * @brief Compute doublet AIC and source AIC from body panel ej to body panel ei
+ */
+std::vector<sdpmDouble> Builder::computeAIC(bool symy, Element const *ei, Element const *ej, Eigen::MatrixXd const &sCoords, sdpmMatrix3d const &g2l)
+{
+ // Compute distance between panels and target panel coordinates in local axes
+ sdpmVector3d cgSrc = ej->getCg();
+ sdpmVector3d cgTgt = ei->getCg();
+ sdpmVector3d distance = cgTgt - cgSrc;
+ if (symy)
+ distance(1) = cgTgt(1) + cgSrc(1);
+ distance = g2l * distance;
+ sdpmVector3d tCoords = g2l * ei->getCg();
+
+ // Compute coefficients
+ int nV = static_cast<int>(ej->getNodes().size());
+ Eigen::Vector4d d, r;
+ for (int i = 1; i <= nV; ++i)
+ {
+ int j = i % nV + 1;
+ d(i - 1) = sqrt((sCoords(0, j - 1) - sCoords(0, i - 1)) * (sCoords(0, j - 1) - sCoords(0, i - 1)) + (sCoords(1, j - 1) - sCoords(1, i - 1)) * (sCoords(1, j - 1) - sCoords(1, i - 1)));
+ r(i - 1) = sqrt((tCoords(0) - sCoords(0, i - 1)) * (tCoords(0) - sCoords(0, i - 1)) + (tCoords(1) - sCoords(1, i - 1)) * (tCoords(1) - sCoords(1, i - 1)) + distance(2) * distance(2));
+ }
+
+ // Compute doublet AIC
+ sdpmDouble mu = 0;
+ if (ei->getId() == ej->getId() && !symy)
+ mu = 0.5; // self-influence
+ else
+ {
+ // additional coefficients
+ Eigen::Vector4d e, h, F, G;
+ for (int i = 1; i <= nV; ++i)
+ {
+ e(i - 1) = (tCoords(0) - sCoords(0, i - 1)) * (tCoords(0) - sCoords(0, i - 1)) + distance(2) * distance(2);
+ h(i - 1) = (tCoords(0) - sCoords(0, i - 1)) * (tCoords(1) - sCoords(1, i - 1));
+ }
+ for (int i = 1; i <= nV; ++i)
+ {
+ int j = i % nV + 1;
+ F(i - 1) = (sCoords(1, j - 1) - sCoords(1, i - 1)) * e(i - 1) - (sCoords(0, j - 1) - sCoords(0, i - 1)) * h(i - 1);
+ G(i - 1) = (sCoords(1, j - 1) - sCoords(1, i - 1)) * e(j - 1) - (sCoords(0, j - 1) - sCoords(0, i - 1)) * h(j - 1);
+ }
+ // AIC
+ for (int i = 1; i <= nV; ++i)
+ {
+ int j = i % nV + 1;
+ mu += atan2(distance(2) * (sCoords(0, j - 1) - sCoords(0, i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)),
+ distance(2) * distance(2) * (sCoords(0, j - 1) - sCoords(0, i - 1)) * (sCoords(0, j - 1) - sCoords(0, i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1));
+ }
+ mu *= -1 / (4 * 3.14159);
+ }
+
+ // Compute source AIC
+ sdpmDouble tau = 0;
+ if (ei->getId() == ej->getId() && !symy)
+ {
+ for (int i = 1; i <= nV; ++i)
+ {
+ int j = i % nV + 1;
+ tau += ((tCoords(0) - sCoords(0, i - 1)) * (sCoords(1, j - 1) - sCoords(1, i - 1)) - (tCoords(1) - sCoords(1, i - 1)) * (sCoords(0, j - 1) - sCoords(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 * 3.14159); // self-influence
+ }
+ else
+ {
+ for (int i = 1; i <= nV; ++i)
+ {
+ int j = i % nV + 1;
+ tau += ((tCoords(0) - sCoords(0, i - 1)) * (sCoords(1, j - 1) - sCoords(1, i - 1)) - (tCoords(1) - sCoords(1, i - 1)) * (sCoords(0, j - 1) - sCoords(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 * 3.14159) * tau - distance(2) * mu;
+ }
+ return std::vector<sdpmDouble>{mu, tau};
+}
diff --git a/sdpm/src/sdpmBuilder.h b/sdpm/src/sdpmBuilder.h
new file mode 100644
index 0000000..d569343
--- /dev/null
+++ b/sdpm/src/sdpmBuilder.h
@@ -0,0 +1,40 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMBUILDER_H
+#define SDPMBUILDER_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include <vector>
+
+namespace sdpm
+{
+
+/**
+ * @brief Aerodynamic Influence Coefficients matrices builder
+ * @authors Adrien Crovato
+ */
+class SDPM_API Builder
+{
+public:
+ static sdpmMatrix3d computeTransfoMat(bool symy, Element const *ej);
+ static sdpmMatrixXd gatherCoords(bool symy, Element const *ej);
+ static std::vector<sdpmDouble> computeAIC(bool symy, Element const *ei, Element const *ej, sdpmMatrixXd const &trsfCoordSrc, sdpmMatrix3d const &g2l);
+};
+
+} // namespace sdpm
+#endif // SDPMBUILDER_H
diff --git a/sdpm/src/sdpmEdge.h b/sdpm/src/sdpmEdge.h
new file mode 100644
index 0000000..41450d8
--- /dev/null
+++ b/sdpm/src/sdpmEdge.h
@@ -0,0 +1,91 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMEDGE_H
+#define SDPMEDGE_H
+
+#include "sdpm.h"
+#include "sdpmNode.h"
+
+#ifndef SWIG
+namespace sdpm
+{
+
+/**
+ * @brief Common edge of three elements
+ * @authors Adrien Crovato
+ */
+class SDPM_API Edge
+{
+public:
+ std::vector<Node *> nods;
+ Element *el0;
+ Element *el1;
+ Element *el2;
+
+ Edge(std::vector<Node *> const &_nods) : nods(_nods), el0(NULL), el1(NULL), el2(NULL) {}
+};
+
+/**
+ * @brief Edge comparator
+ * @authors Adrien Crovato
+ */
+class SDPM_API EquEdge
+{
+public:
+ // Default parameters
+ static const size_t bucket_size = 4;
+ static const size_t min_buckets = 8;
+ // Function returning a unique size_t per Edge
+ size_t operator()(Edge *const f) const
+ {
+ size_t sum = 0;
+ for (auto n : f->nods)
+ sum += static_cast<size_t>(n->getId());
+
+ return sum; // sum of nodes id
+ }
+ // Function allowing to univoquely identify a Edge
+ bool operator()(Edge *const f0, Edge *const f1) const
+ {
+ bool flag = false;
+ size_t cnt = 0;
+ // compare nodes of f0 to nodes of f1
+ for (auto n0 : f0->nods)
+ {
+ for (auto n1 : f1->nods)
+ {
+ if (n0 == n1)
+ {
+ cnt++;
+ break;
+ }
+ }
+ if (cnt == 0)
+ break;
+ }
+ // check the number of shared nodes
+ if (cnt == f0->nods.size())
+ flag = true;
+
+ return flag; // true if f0 = f1 (ie, f0 and f1 share the same nodes)
+ }
+};
+
+} // namespace sdpm
+#endif
+
+#endif // SDPMEDGE_H
diff --git a/sdpm/src/sdpmElement.cpp b/sdpm/src/sdpmElement.cpp
index a795dd6..53b092f 100644
--- a/sdpm/src/sdpmElement.cpp
+++ b/sdpm/src/sdpmElement.cpp
@@ -20,7 +20,6 @@ using namespace sdpm;
namespace sdpm
{
-
SDPM_API std::ostream &operator<<(std::ostream &out, ElType const &obj)
{
switch (obj)
@@ -37,41 +36,28 @@ SDPM_API std::ostream &operator<<(std::ostream &out, ElType const &obj)
}
return out;
}
-
} // namespace sdpm
-Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes), _hasVals(false)
+Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes)
{
_tag->addElements(this);
}
-void Element::init()
-{
- throw std::runtime_error("Element::init() not implemented!\n");
-}
-void Element::update()
-{
- throw std::runtime_error("Element::update() not implemented!\n");
-}
-
/**
- * @brief Compute center of gravity of element
+ * @brief Update precomputed values
*/
-void Element::computeCg()
+void Element::update()
{
+ // Center of gravity
_cg = Eigen::Vector3d::Zero();
for (auto n : _nodes)
- _cg += n->getCoord();
+ _cg += n->getCoords();
_cg /= static_cast<double>(_nodes.size());
}
-void Element::computeVolume()
-{
- throw std::runtime_error("Element::computeVolume not implemented\n");
-}
-void Element::computeNormal()
+sdpmVector3d Element::computeGradient(std::vector<sdpmDouble> const &u) const
{
- throw std::runtime_error("Element::computeNormal not implemented\n");
+ throw std::runtime_error("Element::computeGradient not implemented\n");
}
void Element::write(std::ostream &out) const
diff --git a/sdpm/src/sdpmElement.h b/sdpm/src/sdpmElement.h
index 0a32e32..938aad7 100644
--- a/sdpm/src/sdpmElement.h
+++ b/sdpm/src/sdpmElement.h
@@ -51,8 +51,6 @@ protected:
size_t _id; ///< element number
Tag *_tag; ///< physical tag to which element belongs
std::vector<Node *> _nodes; ///< nodes of the element
- // Flags
- bool _hasVals; ///< true if (precomputed) values have been initialized
// Geometry
sdpmVector3d _cg; ///< element center of gravity
sdpmDouble _vol; ///< element volume
@@ -70,32 +68,24 @@ public:
}
#ifndef SWIG
-protected:
- // Geometry
- void computeCg();
- virtual void computeVolume();
- virtual void computeNormal();
-
-public:
// Management
- virtual void init();
virtual void update();
+ // Setters
+ void setNode(size_t i, Node *n) { _nodes[i] = n; }
// Getters
size_t getId() const { return _id; }
size_t const &getIdRef() const { return _id; }
Tag const *getTag() const { return _tag; }
std::vector<Node *> const &getNodes() const { return _nodes; }
- inline bool hasValues() const;
- inline sdpmVector3d const &getCg() const;
- inline sdpmDouble getVolume() const;
- inline sdpmVector3d const &getNormal() const;
+ sdpmVector3d const &getCg() const { return _cg; }
+ sdpmDouble getVolume() const { return _vol; }
+ sdpmVector3d const &getNormal() const { return _nrm; }
// Utilities
+ virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const;
virtual void write(std::ostream &out) const override;
#endif
};
-#include "sdpmElement.inl"
-
} // namespace sdpm
#endif // SDPMELEMENT_H
diff --git a/sdpm/src/sdpmElement.inl b/sdpm/src/sdpmElement.inl
deleted file mode 100644
index b1b74f9..0000000
--- a/sdpm/src/sdpmElement.inl
+++ /dev/null
@@ -1,56 +0,0 @@
-/*
- * Copyright 2023 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.
- */
-
-/**
- * @brief Return wether the element values have been initialized or not
- */
-inline bool Element::hasValues() const
-{
- return _hasVals;
-}
-
-/**
- * @brief Return the volume of the element
- */
-inline sdpmDouble Element::getVolume() const
-{
- if (_hasVals)
- return _vol;
- else
- throw std::runtime_error("Element::getVolume is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the center of gravity of the element
- */
-inline sdpmVector3d const &Element::getCg() const
-{
- if (_hasVals)
- return _cg;
- else
- throw std::runtime_error("Element::getCg is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the unit normal vector of the element
- */
-inline sdpmVector3d const &Element::getNormal() const
-{
- if (_hasVals)
- return _nrm;
- else
- throw std::runtime_error("Element::getNormal is invalid because element memory is not up-to-date!\n");
-}
diff --git a/sdpm/src/sdpmGmshExport.cpp b/sdpm/src/sdpmGmshExport.cpp
index 987c274..c83af97 100644
--- a/sdpm/src/sdpmGmshExport.cpp
+++ b/sdpm/src/sdpmGmshExport.cpp
@@ -29,10 +29,10 @@ GmshExport::~GmshExport() { std::cout << "~sdpm::GmshExport()\n"; }
/**
* @brief Save the mesh on disk
*/
-void GmshExport::save(std::string const &fname) const
+void GmshExport::save(std::string const &sfx) const
{
// write gmsh .msh file with stdio commands
- std::string wname = fname + ".msh";
+ std::string wname = _msh.getName() + sfx + ".msh";
std::cout << "writing file: " << wname << "... " << std::flush;
FILE *file = fopen(wname.c_str(), "wt");
if (!file)
@@ -44,26 +44,29 @@ void GmshExport::save(std::string const &fname) const
fprintf(file, "$EndMeshFormat\n");
// write physical names
+ std::map<std::string, Tag *> tags = _msh.getTags();
fprintf(file, "$PhysicalNames\n");
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.tags.size()));
- for (auto &p : _msh.tags)
+ fprintf(file, "%zu\n", static_cast<size_t>(tags.size()));
+ for (auto &p : tags)
fprintf(file, "%zu %zu \"%s\"\n", p.second->getDim(), p.second->getId(), p.second->getName().c_str());
fprintf(file, "$EndPhysicalNames\n");
// nodes
+ std::vector<Node *> nods = _msh.getNodes();
fprintf(file, "$Nodes\n");
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.nodes.size()));
- for (auto n : _msh.nodes)
+ fprintf(file, "%zu\n", static_cast<size_t>(nods.size()));
+ for (auto n : nods)
{
- sdpmVector3d pos = n->getCoord();
+ sdpmVector3d pos = n->getCoords();
fprintf(file, "%zu %.16g %.16g %.16g\n", n->getId(), pos(0), pos(1), pos(2));
}
fprintf(file, "$EndNodes\n");
// elements
+ std::vector<Element *> elems = _msh.getElements();
fprintf(file, "$Elements\n");
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.elements.size()));
- for (auto e : _msh.elements)
+ fprintf(file, "%zu\n", static_cast<size_t>(elems.size()));
+ for (auto e : elems)
{
fprintf(file, "%zu %zu %zu", e->getId(), static_cast<size_t>(e->getType()), static_cast<size_t>(2));
// physical tag
@@ -82,16 +85,16 @@ void GmshExport::save(std::string const &fname) const
// close file
fclose(file);
- std::cout << "done" << std::endl;
+ std::cout << "done." << std::endl;
}
/**
* @brief Save results (associated to a mesh) on disk
*/
-void GmshExport::save(std::string const &fname, Results const &r) const
+void GmshExport::save(Results const &r, std::string const &sfx) const
{
// write results to file
- std::string wname = fname + ".pos";
+ std::string wname = _msh.getName() + sfx + ".pos";
std::cout << "writing file: " << wname << "... " << std::flush;
FILE *file;
if (_bin)
@@ -111,27 +114,31 @@ void GmshExport::save(std::string const &fname, Results const &r) const
}
fprintf(file, "$EndMeshFormat\n");
+ // get nodes and elements
+ std::vector<Node *> nods = _msh.getNodes();
+ std::vector<Element *> elems = _msh.getElements();
+
// node data (scalar)
for (auto &p : r.scalarsAtNodes)
{
fprintf(file, "$NodeData\n");
- fprintf(file, "%d\n", 1); // nb of string tags
- fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
- fprintf(file, "%d\n", 1); // nb of real tags (1)
- fprintf(file, "%.16g\n", r.time); // 1. the time value
- fprintf(file, "%d\n", 3); // nb of integer tags (3)
- fprintf(file, "%d\n", r.nt); // 1. time step index
- fprintf(file, "%d\n", 1); // 2. nb of field component (1, 3 or 9)
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.nodes.size())); // 3. nb of entities in the view
- for (size_t ii = 0; ii < _msh.nodes.size(); ++ii)
+ fprintf(file, "%d\n", 1); // nb of string tags
+ fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
+ fprintf(file, "%d\n", 1); // nb of real tags (1)
+ fprintf(file, "%.16g\n", r.time); // 1. the time value
+ fprintf(file, "%d\n", 3); // nb of integer tags (3)
+ fprintf(file, "%d\n", r.nt); // 1. time step index
+ fprintf(file, "%d\n", 1); // 2. nb of field component (1, 3 or 9)
+ fprintf(file, "%zu\n", static_cast<size_t>(nods.size())); // 3. nb of entities in the view
+ for (size_t ii = 0; ii < nods.size(); ++ii)
{
if (_bin)
{
- fwrite(&(_msh.nodes[ii]->getIdRef()), sizeof(int), 1, file);
+ fwrite(&(nods[ii]->getIdRef()), sizeof(int), 1, file);
fwrite(&((*(p.second))[ii]), sizeof(sdpmDouble), 1, file);
}
else
- fprintf(file, "%zu %.16g\n", _msh.nodes[ii]->getId(), (*(p.second))[ii]);
+ fprintf(file, "%zu %.16g\n", nods[ii]->getId(), (*(p.second))[ii]);
}
fprintf(file, "$EndNodeData\n");
}
@@ -140,23 +147,23 @@ void GmshExport::save(std::string const &fname, Results const &r) const
for (auto &p : r.scalarsAtElems)
{
fprintf(file, "$ElementData\n");
- fprintf(file, "%d\n", 1); // nb of string tags
- fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
- fprintf(file, "%d\n", 1); // nb of real tags (1)
- fprintf(file, "%.16g\n", r.time); // 1. the time value
- fprintf(file, "%d\n", 3); // nb of integer tags (3)
- fprintf(file, "%d\n", r.nt); // 1. time step index
- fprintf(file, "%d\n", 1); // 2. nb of field component (1, 3 or 9)
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.elements.size())); // 3. nb of entities in the view
- for (size_t ii = 0; ii < _msh.elements.size(); ++ii)
+ fprintf(file, "%d\n", 1); // nb of string tags
+ fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
+ fprintf(file, "%d\n", 1); // nb of real tags (1)
+ fprintf(file, "%.16g\n", r.time); // 1. the time value
+ fprintf(file, "%d\n", 3); // nb of integer tags (3)
+ fprintf(file, "%d\n", r.nt); // 1. time step index
+ fprintf(file, "%d\n", 1); // 2. nb of field component (1, 3 or 9)
+ fprintf(file, "%zu\n", static_cast<size_t>(elems.size())); // 3. nb of entities in the view
+ for (size_t ii = 0; ii < elems.size(); ++ii)
{
if (_bin)
{
- fwrite(&(_msh.elements[ii]->getIdRef()), sizeof(int), 1, file);
+ fwrite(&(elems[ii]->getIdRef()), sizeof(int), 1, file);
fwrite(&(*(p.second))[ii], sizeof(sdpmDouble), 1, file);
}
else
- fprintf(file, "%zu %.16g\n", _msh.elements[ii]->getId(), (*(p.second))[ii]);
+ fprintf(file, "%zu %.16g\n", elems[ii]->getId(), (*(p.second))[ii]);
}
fprintf(file, "$EndElementData\n");
}
@@ -165,24 +172,24 @@ void GmshExport::save(std::string const &fname, Results const &r) const
for (auto &p : r.vectorsAtNodes)
{
fprintf(file, "$NodeData\n");
- fprintf(file, "%d\n", 1); // nb of string tags
- fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
- fprintf(file, "%d\n", 1); // nb of real tags (1)
- fprintf(file, "%.16g\n", r.time); // 1. the time value
- fprintf(file, "%d\n", 3); // nb of integer tags (3)
- fprintf(file, "%d\n", r.nt); // 1. time step index
- fprintf(file, "%d\n", 3); // 2. nb of field component (1, 3 or 9)
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.nodes.size())); // 3. nb of entities in the view
- for (size_t ii = 0; ii < _msh.nodes.size(); ++ii)
+ fprintf(file, "%d\n", 1); // nb of string tags
+ fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
+ fprintf(file, "%d\n", 1); // nb of real tags (1)
+ fprintf(file, "%.16g\n", r.time); // 1. the time value
+ fprintf(file, "%d\n", 3); // nb of integer tags (3)
+ fprintf(file, "%d\n", r.nt); // 1. time step index
+ fprintf(file, "%d\n", 3); // 2. nb of field component (1, 3 or 9)
+ fprintf(file, "%zu\n", static_cast<size_t>(nods.size())); // 3. nb of entities in the view
+ for (size_t ii = 0; ii < nods.size(); ++ii)
{
sdpmVector3d const &v = (*(p.second))[ii];
if (_bin)
{
- fwrite(&(_msh.nodes[ii]->getIdRef()), sizeof(int), 1, file);
+ fwrite(&(nods[ii]->getIdRef()), sizeof(int), 1, file);
fwrite(&(v(0)), sizeof(sdpmDouble), 3, file);
}
else
- fprintf(file, "%zu %.16g %.16g %.16g\n", _msh.nodes[ii]->getId(), v(0), v(1), v(2));
+ fprintf(file, "%zu %.16g %.16g %.16g\n", nods[ii]->getId(), v(0), v(1), v(2));
}
fprintf(file, "$EndNodeData\n");
}
@@ -191,24 +198,24 @@ void GmshExport::save(std::string const &fname, Results const &r) const
for (auto &p : r.vectorsAtElems)
{
fprintf(file, "$ElementData\n");
- fprintf(file, "%d\n", 1); // nb of string tags
- fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
- fprintf(file, "%d\n", 1); // nb of real tags (1)
- fprintf(file, "%.16g\n", r.time); // 1. the time value
- fprintf(file, "%d\n", 3); // nb of integer tags (3)
- fprintf(file, "%d\n", r.nt); // 1. time step index
- fprintf(file, "%d\n", 3); // 2. nb of field component (1, 3 or 9)
- fprintf(file, "%zu\n", static_cast<size_t>(_msh.elements.size())); // 3. nb of entities in the view
- for (size_t ii = 0; ii < _msh.elements.size(); ++ii)
+ fprintf(file, "%d\n", 1); // nb of string tags
+ fprintf(file, "\"%s\"\n", p.first.c_str()); // the name of the view
+ fprintf(file, "%d\n", 1); // nb of real tags (1)
+ fprintf(file, "%.16g\n", r.time); // 1. the time value
+ fprintf(file, "%d\n", 3); // nb of integer tags (3)
+ fprintf(file, "%d\n", r.nt); // 1. time step index
+ fprintf(file, "%d\n", 3); // 2. nb of field component (1, 3 or 9)
+ fprintf(file, "%zu\n", static_cast<size_t>(elems.size())); // 3. nb of entities in the view
+ for (size_t ii = 0; ii < elems.size(); ++ii)
{
sdpmVector3d const &v = (*(p.second))[ii];
if (_bin)
{
- fwrite(&(_msh.elements[ii]->getIdRef()), sizeof(int), 1, file);
+ fwrite(&(elems[ii]->getIdRef()), sizeof(int), 1, file);
fwrite(&(v(0)), sizeof(sdpmDouble), 3, file);
}
else
- fprintf(file, "%zu %.16g %.16g %.16g\n", _msh.elements[ii]->getId(), v(0), v(1), v(2));
+ fprintf(file, "%zu %.16g %.16g %.16g\n", elems[ii]->getId(), v(0), v(1), v(2));
}
fprintf(file, "$EndElementData\n");
}
diff --git a/sdpm/src/sdpmGmshExport.h b/sdpm/src/sdpmGmshExport.h
index 8e9051f..6ecbabf 100644
--- a/sdpm/src/sdpmGmshExport.h
+++ b/sdpm/src/sdpmGmshExport.h
@@ -30,6 +30,7 @@ namespace sdpm
/**
* @brief Write Gmsh mesh and result files
* @authors Adrien Crovato
+ * @todo replace by Gmsh API?
*/
class SDPM_API GmshExport : public sdpmObject
{
@@ -37,11 +38,11 @@ class SDPM_API GmshExport : public sdpmObject
bool _bin; ///< binary flag
public:
GmshExport(Mesh &msh, bool b = true);
- virtual ~GmshExport();
+ ~GmshExport();
- void save(std::string const &fname) const;
+ void save(std::string const &sfx = "") const;
#ifndef SWIG
- void save(std::string const &fname, Results const &r) const;
+ void save(Results const &r, std::string const &sfx = "") const;
#endif
virtual void write(std::ostream &out) const override;
};
diff --git a/sdpm/src/sdpmGmshImport.cpp b/sdpm/src/sdpmGmshImport.cpp
index 2d6bef0..7cfd10f 100644
--- a/sdpm/src/sdpmGmshImport.cpp
+++ b/sdpm/src/sdpmGmshImport.cpp
@@ -84,12 +84,13 @@ void GmshImport::readPhysicalTags(FILE *file)
if (fscanf(file, "%zu %zu \"%[^\"]\"", &dim, &no, name) != 3)
throw std::runtime_error("GmshImport::readPhysicalTags: error reading dim, no, name!\n");
// add tag to tag list
- auto it = _msh.tags.find(name);
- if (it == _msh.tags.end())
+ std::map<std::string, Tag *> tags = _msh.getTags();
+ auto it = tags.find(name);
+ if (it == tags.end())
{
Tag *tag = new Tag(no, name, dim);
_tagMap[no] = tag;
- _msh.tags[name] = tag;
+ _msh.addTag(static_cast<std::string>(name), tag);
}
}
}
@@ -106,7 +107,7 @@ void GmshImport::readNodes(FILE *file)
std::cout << "\tcounted " << nbnod << " nodes\n";
// allocate memory and read the nodes
- _msh.nodes.reserve(nbnod);
+ _msh.reserveNodes(nbnod);
for (int i = 0; i < nbnod; ++i)
{
size_t no = 0;
@@ -116,8 +117,7 @@ void GmshImport::readNodes(FILE *file)
if (_verb)
std::cout << "\t" << no << ((i == nbnod - 1) ? ".\n" : ",");
Node *nod = new Node(no, sdpmVector3d(x, y, z));
- _msh.nodes.push_back(nod);
-
+ _msh.addNode(nod);
_nodeMap[no] = nod;
}
}
@@ -134,7 +134,7 @@ void GmshImport::readElements(FILE *file)
std::cout << "\tcounted " << nbelm << " elements\n";
// allocate memory and read the elements
- _msh.elements.reserve(nbelm);
+ _msh.reserveElements(nbelm);
for (int i = 0; i < nbelm; ++i)
{
// read number, type and nb of tags
@@ -161,12 +161,12 @@ void GmshImport::readElements(FILE *file)
if (static_cast<ElType>(type) == ElType::QUAD4) // 4-node quad
{
readElNods(file, 4, nods);
- _msh.elements.push_back(new Quad4(no, ptag, nods));
+ _msh.addElement(new Quad4(no, ptag, nods));
}
else if (static_cast<ElType>(type) == ElType::LINE2) // 2-node line
{
readElNods(file, 2, nods);
- _msh.elements.push_back(new Line2(no, ptag, nods));
+ _msh.addElement(new Line2(no, ptag, nods));
}
else
skipline = true;
@@ -234,8 +234,8 @@ void GmshImport::load(std::string const &fname)
auto idx = fname.find(".msh");
if (idx != std::string::npos)
{
- _msh.name = fname.substr(0, idx);
- std::cout << "\tsetting mesh name to \"" << _msh.name << "\"\n";
+ _msh.setName(fname.substr(0, idx));
+ std::cout << "\tsetting mesh name to \"" << _msh.getName() << "\"\n";
}
// read gmsh .msh file with stdio commands
diff --git a/sdpm/src/sdpmGmshImport.h b/sdpm/src/sdpmGmshImport.h
index 8a9466f..9d07ad4 100644
--- a/sdpm/src/sdpmGmshImport.h
+++ b/sdpm/src/sdpmGmshImport.h
@@ -29,7 +29,7 @@ namespace sdpm
/**
* @brief Read Gmsh mesh file and load into data structure
* @authors Adrien Crovato
- * @todo replace by Gmsh API
+ * @todo replace by Gmsh API?
*/
class SDPM_API GmshImport : public sdpmObject
{
@@ -48,7 +48,7 @@ class SDPM_API GmshImport : public sdpmObject
public:
GmshImport(Mesh &msh, bool v = false);
- virtual ~GmshImport();
+ ~GmshImport();
void load(std::string const &fname);
void write(std::ostream &out) const;
diff --git a/sdpm/src/sdpmGroup.cpp b/sdpm/src/sdpmGroup.cpp
new file mode 100644
index 0000000..b61848c
--- /dev/null
+++ b/sdpm/src/sdpmGroup.cpp
@@ -0,0 +1,40 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpmGroup.h"
+#include "sdpmMesh.h"
+
+using namespace sdpm;
+
+Group::Group(Mesh &msh, std::string const &name) : sdpmObject(), _msh(msh)
+{
+ // Check if the group named "name" exists
+ std::map<std::string, Tag *> tags = _msh.getTags();
+ auto it = tags.find(name);
+ if (it == tags.end())
+ {
+ std::stringstream out;
+ out << "sdpm::Group: Physical tag \"" << name << "\" not found among available tags:\n";
+ for (auto it = tags.begin(); it != tags.end(); ++it)
+ {
+ auto itn(it);
+ ++itn;
+ out << (it->first) << ((itn != tags.end()) ? ", " : "\n\n");
+ }
+ throw std::out_of_range(out.str());
+ }
+ _tag = it->second;
+}
diff --git a/sdpm/src/sdpmGroup.h b/sdpm/src/sdpmGroup.h
new file mode 100644
index 0000000..ce859bd
--- /dev/null
+++ b/sdpm/src/sdpmGroup.h
@@ -0,0 +1,48 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMGROUP_H
+#define SDPMGROUP_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include "sdpmTag.h"
+#include <vector>
+#include <string>
+
+namespace sdpm
+{
+
+/**
+ * @brief Group of elements
+ * @authors Adrien Crovato
+ */
+class SDPM_API Group : public sdpmObject
+{
+protected:
+ Mesh &_msh; ///< mesh data structure
+ Tag *_tag; ///< named (physical) tag holding the elements
+
+public:
+ Group(Mesh &msh, std::string const &name);
+
+ std::vector<Element *> const &getElements() const { return _tag->getElements(); }
+ std::string const &getName() const { return _tag->getName(); }
+};
+
+} // namespace sdpm
+
+#endif // SDPMGROUP_H
diff --git a/sdpm/src/sdpmLine2.cpp b/sdpm/src/sdpmLine2.cpp
index d20bc84..21f7b6b 100644
--- a/sdpm/src/sdpmLine2.cpp
+++ b/sdpm/src/sdpmLine2.cpp
@@ -18,14 +18,8 @@
using namespace sdpm;
-/**
- * @brief Initialize precomputed values
- */
-void Line2::init()
+Line2::Line2(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
{
- // Set flags
- _hasVals = true;
- // Update values
update();
}
@@ -34,29 +28,7 @@ void Line2::init()
*/
void Line2::update()
{
- // Update values
- if (_hasVals)
- {
- computeCg();
- //computeVolume();
- //computeNormal();
- }
-}
-
-/**
- * @brief Return the element surface
- */
-void Line2::computeVolume()
-{
- throw std::runtime_error("Not Implemented yet.\n");
- // TODO
-}
-
-/**
- * @brief Compute element unit normal vector
- */
-void Line2::computeNormal()
-{
- throw std::runtime_error("Not Implemented yet.\n");
- // TODO
+ Element::update(); // cg
+ _vol = (_nodes[1]->getCoords() - _nodes[0]->getCoords()).norm(); // volume
+ _nrm = sdpmVector3d::Zero(); // unit normal
}
diff --git a/sdpm/src/sdpmLine2.h b/sdpm/src/sdpmLine2.h
index 014287b..8eea05e 100644
--- a/sdpm/src/sdpmLine2.h
+++ b/sdpm/src/sdpmLine2.h
@@ -29,19 +29,12 @@ namespace sdpm
*/
class SDPM_API Line2 : public Element
{
-#ifndef SWIG
-protected:
- virtual void computeVolume() override;
- virtual void computeNormal() override;
-#endif
-
public:
- Line2(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes) {}
- virtual ~Line2() {}
+ Line2(size_t n, Tag *tag, std::vector<Node *> &nodes);
+ ~Line2() {}
virtual ElType getType() const override { return ElType::LINE2; }
#ifndef SWIG
- virtual void init() override;
virtual void update() override;
#endif
};
diff --git a/sdpm/src/sdpmMesh.cpp b/sdpm/src/sdpmMesh.cpp
index 41c4fe2..5c2828f 100644
--- a/sdpm/src/sdpmMesh.cpp
+++ b/sdpm/src/sdpmMesh.cpp
@@ -22,7 +22,7 @@
using namespace sdpm;
-Mesh::Mesh() : sdpmObject(), name("none") {}
+Mesh::Mesh() : sdpmObject(), _name("none") {}
Mesh::~Mesh()
{
@@ -35,21 +35,21 @@ Mesh::~Mesh()
*/
void Mesh::clear()
{
- for (auto n : nodes)
+ for (auto n : _nodes)
delete n;
- nodes.clear();
- for (auto e : elements)
+ _nodes.clear();
+ for (auto e : _elements)
delete e;
- elements.clear();
- for (auto &t : tags)
+ _elements.clear();
+ for (auto &t : _tags)
delete t.second;
- tags.clear();
+ _tags.clear();
}
void Mesh::write(std::ostream &out) const
{
- out << "sdpmMesh \"" << name << "\": \n";
- out << "\t" << nodes.size() << " nodes\n";
- out << "\t" << elements.size() << " elements\n";
- out << "\t" << tags.size() << " physical tags";
+ out << "sdpmMesh \"" << _name << "\": \n";
+ out << "\t" << _nodes.size() << " nodes\n";
+ out << "\t" << _elements.size() << " elements\n";
+ out << "\t" << _tags.size() << " physical tags";
}
diff --git a/sdpm/src/sdpmMesh.h b/sdpm/src/sdpmMesh.h
index adae178..a85aa8a 100644
--- a/sdpm/src/sdpmMesh.h
+++ b/sdpm/src/sdpmMesh.h
@@ -35,14 +35,28 @@ namespace sdpm
*/
class SDPM_API Mesh : public sdpmObject
{
+ std::string _name; ///< mesh name
+ std::vector<Node *> _nodes; ///< nodes
+ std::vector<Element *> _elements; ///< elements
+ std::map<std::string, Tag *> _tags; ///< named (physical) tags
public:
- std::string name; ///< mesh name
- std::vector<Node *> nodes; ///< nodes
- std::vector<Element *> elements; ///< elements
- std::map<std::string, Tag *> tags; ///< named (physical) tags
-
Mesh();
- virtual ~Mesh();
+ ~Mesh();
+
+ // Getters
+ std::string const &getName() const { return _name; }
+ std::vector<Node *> const &getNodes() const { return _nodes; }
+ std::vector<Element *> const &getElements() const { return _elements; }
+ std::map<std::string, Tag *> const &getTags() const { return _tags; }
+#ifndef SWIG
+ // Setters
+ void setName(std::string &n) { _name = n; }
+ void addNode(Node *n) { _nodes.push_back(n); }
+ void addElement(Element *e) { _elements.push_back(e); }
+ void addTag(std::string &n, Tag *t) { _tags[n] = t; }
+ void reserveNodes(size_t n) { _nodes.reserve(n); }
+ void reserveElements(size_t n) { _elements.reserve(n); }
+#endif
void clear();
#ifndef SWIG
diff --git a/sdpm/src/sdpmNode.cpp b/sdpm/src/sdpmNode.cpp
index c62f666..3d4f5e4 100644
--- a/sdpm/src/sdpmNode.cpp
+++ b/sdpm/src/sdpmNode.cpp
@@ -15,10 +15,9 @@
*/
#include "sdpmNode.h"
-
using namespace sdpm;
void Node::write(std::ostream &out) const
{
- out << "node #" << _id << " = " << _coord.transpose();
+ out << "node #" << _id << " = [" << _coord.transpose() << "]";
}
diff --git a/sdpm/src/sdpmNode.h b/sdpm/src/sdpmNode.h
index 3a242b9..5fc106b 100644
--- a/sdpm/src/sdpmNode.h
+++ b/sdpm/src/sdpmNode.h
@@ -37,7 +37,7 @@ public:
size_t getId() { return _id; }
size_t const &getIdRef() { return _id; }
- sdpmVector3d getCoord() { return _coord; }
+ sdpmVector3d getCoords() { return _coord; }
#ifndef SWIG
virtual void write(std::ostream &out) const override;
diff --git a/sdpm/src/sdpmProblem.cpp b/sdpm/src/sdpmProblem.cpp
new file mode 100644
index 0000000..ad141a3
--- /dev/null
+++ b/sdpm/src/sdpmProblem.cpp
@@ -0,0 +1,112 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include "sdpmElement.h"
+#include "sdpmTag.h"
+#include <iostream>
+using namespace sdpm;
+
+Problem::Problem(Mesh &msh, double aoa, double aos, double mch,
+ double sref, double cref, double xref, double yref, double zref,
+ bool ysym) : _msh(msh), _mach(mch), _sRef(sref), _cRef(cref), _ySym(ysym)
+{
+ updateFreestream(aoa, aos);
+ _xRef(0) = xref;
+ _xRef(1) = yref;
+ _xRef(2) = zref;
+}
+
+Problem::~Problem()
+{
+ std::cout << "~sdpm::Problem()\n";
+}
+
+/**
+ * @brief Add a non-lifting body
+ */
+void Problem::addBody(Body &b)
+{
+ _bodies.push_back(&b);
+}
+
+/**
+ * @brief Update the freestream values
+ */
+void Problem::updateFreestream(double aoa, double aos)
+{
+ // Angles and velocity
+ _aoa = aoa;
+ _aos = aos;
+ _uinf = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
+ // Directions
+ _dDir = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
+ _sDir = sdpmVector3d(-cos(_aoa) * sin(_aos), cos(_aos), -sin(_aoa) * sin(_aos));
+ _lDir = sdpmVector3d(-sin(_aoa), 0, cos(_aoa));
+}
+
+/**
+ * @brief Update the elements values
+ */
+void Problem::updateElements()
+{
+ // Update surface Jacobian
+ for (auto b : _bodies)
+ {
+ for (auto e : b->getElements())
+ e->update();
+ for (auto e : b->getWake().getElements())
+ e->update();
+ }
+}
+
+/**
+ * @brief Check that Problem is not empty and that elements are supported
+ */
+void Problem::check() const
+{
+ // Sanity checks
+ if (_bodies.empty())
+ throw std::runtime_error("sdpm::Problem: no lifting bodies provided!\n");
+ if (_mach != 0)
+ throw std::runtime_error("sdpm::Problem: compressible computation feature not implemented yet!\n");
+
+ // Surface elements
+ for (auto b : _bodies)
+ for (auto e : b->getElements())
+ if (e->getType() != ElType::QUAD4)
+ {
+ std::stringstream err;
+ err << "SDPM is only implemented for surface elements of type "
+ << ElType::QUAD4 << " (" << e->getType() << " was given)!\n";
+ throw std::runtime_error(err.str());
+ }
+}
+
+void Problem::write(std::ostream &out) const
+{
+ out << "sdpm::Problem definition"
+ << "\n\tAngle of attack: " << _aoa
+ << "\n\tAngle of sideslip: " << _aos
+ << "\n\tReference area: " << _sRef
+ << "\n\tReference length: " << _cRef
+ << "\n\tReference point: [" << _xRef.transpose() << "]"
+ << std::endl;
+ out << "with";
+ for (auto b : _bodies)
+ out << "\t" << *b;
+}
diff --git a/sdpm/src/sdpmProblem.h b/sdpm/src/sdpmProblem.h
new file mode 100644
index 0000000..4a2dac8
--- /dev/null
+++ b/sdpm/src/sdpmProblem.h
@@ -0,0 +1,86 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMPROBLEM_H
+#define SDPMPROBLEM_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include "sdpmMesh.h"
+#include <vector>
+
+namespace sdpm
+{
+
+/**
+ * @brief Hold the problem definition
+ * @authors Adrien Crovato
+ */
+class SDPM_API Problem : public sdpmObject
+{
+ // Mesh
+ Mesh &_msh; ///< mesh data structure
+ std::vector<Body *> _bodies; ///< bodies in fluid
+ // Aerodynamic
+ sdpmDouble _aoa; ///< Angle of attack
+ sdpmDouble _aos; ///< Angle of sideslip
+ double _mach; ///< Mach number
+ // Geometry
+ double _sRef; ///< Reference surface area
+ double _cRef; ///< Reference chord length
+ Eigen::Vector3d _xRef; ///< Reference center point for moment computation
+ bool _ySym; ///< True if mesh is symmetric with respect to y-plane
+ // Freestream
+ sdpmVector3d _uinf; ///< freestream velocity
+ sdpmVector3d _dDir; ///< direction of drag force
+ sdpmVector3d _sDir; ///< direction of side force
+ sdpmVector3d _lDir; ///< direction of lift force
+
+public:
+ Problem(Mesh &msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref, bool ysym);
+ ~Problem();
+
+#ifndef SWIG
+ std::vector<Body *> const &getBodies() const
+ {
+ return _bodies;
+ }
+ std::vector<Node *> const &getNodes() const { return _msh.getNodes(); }
+ std::vector<Element *> const &getElements() const { return _msh.getElements(); }
+ sdpmDouble getAoa() const { return _aoa; }
+ sdpmDouble getAos() const { return _aos; }
+ sdpmDouble getRefArea() const { return _sRef; }
+ sdpmDouble getRefLgth() const { return _cRef; }
+ sdpmVector3d const &getRefCtr() const { return _xRef; }
+ sdpmDouble getSym() const { return _ySym; }
+ sdpmVector3d const &getFrstrmVel() const { return _uinf; }
+ sdpmVector3d const &getDragDir() const { return _dDir; }
+ sdpmVector3d const &getSideDir() const { return _sDir; }
+ sdpmVector3d const &getLiftDir() const { return _lDir; }
+#endif
+ void addBody(Body &b);
+ void updateFreestream(double aoa, double aos);
+
+#ifndef SWIG
+ void check() const;
+ void updateElements();
+ virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace sdpm
+
+#endif // SDPMPROBLEM_H
diff --git a/sdpm/src/sdpmQuad4.cpp b/sdpm/src/sdpmQuad4.cpp
index 6caf379..8281830 100644
--- a/sdpm/src/sdpmQuad4.cpp
+++ b/sdpm/src/sdpmQuad4.cpp
@@ -18,49 +18,56 @@
using namespace sdpm;
-/**
- * @brief Initialize precomputed values
- */
-void Quad4::init()
+Quad4::Quad4(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
{
- // Set flags
- _hasVals = true;
- // Update values
update();
}
/**
- * @brief Update precomputed values and gradients
+ * @brief Update precomputed values
*/
void Quad4::update()
{
- // Update values
- if (_hasVals)
- {
- computeCg();
- //computeVolume();
- //computeNormal();
- }
+ Element::update(); // cg
+ _tgt0 = _nodes[2]->getCoords() - _nodes[0]->getCoords(); // unit tangent (n0->n2)
+ _tgt1 = _nodes[3]->getCoords() - _nodes[1]->getCoords(); // unit tangent (n1->n4)
+ _vol = 0.5 * _tgt0.cross(_tgt1).norm(); // volume
+ _nrm = _tgt0.cross(_tgt1).normalized(); // unit normal
}
/**
- * @brief Return the element surface
+ * @brief Compute gradient in element (i.e. surface gradient)
+ * @todo Either precompute dsf once per element type and precompute jacobian once per element, or use finite-volume approach (Gauss theorem)
*/
-void Quad4::computeVolume()
+sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
{
- throw std::runtime_error("Not Implemented yet.\n");
- // TODO use cross products of diagonals
- //std::vector<sdpmVector3d> ts = {(nodes[2]->getCoords() - nodes[0]->getCoords()), (nodes[3]->getCoords() - nodes[1]->getCoords())};
- //return 0.5 * ts[0].cross(ts[1]).norm();
-}
+ // Compute gradient of shape functions
+ sdpmDouble ksi = 0., eta = 0.;
+ sdpmMatrixXd dsf(2, 4);
+ dsf(0, 0) = 0.25 * (1. + eta);
+ dsf(1, 0) = 0.25 * (1. + ksi);
+ dsf(0, 1) = -0.25 * (1. + eta);
+ dsf(1, 1) = 0.25 * (1. - ksi);
+ dsf(0, 2) = -0.25 * (1. - eta);
+ dsf(1, 2) = -0.25 * (1. - ksi);
+ dsf(0, 3) = 0.25 * (1. - eta);
+ dsf(1, 3) = -0.25 * (1. + ksi);
-/**
- * @brief Compute element unit normal vector
- */
-void Quad4::computeNormal()
-{
- throw std::runtime_error("Not Implemented yet.\n");
- // TODO use diagonals or tangents?
- //std::vector<sdpmVector3d> ts = {(nodes[2]->getCoords() - nodes[0]->getCoords()), (nodes[3]->getCoords() - nodes[1]->getCoords())};
- //normal = ts[0].cross(ts[1]).normalized();
+ // Compute the inverse Jacobian in 3D
+ sdpmMatrix3d jac = sdpmMatrix3d::Zero(), ijac = sdpmMatrix3d::Zero();
+ // 1) compute first 2 rows: dXi, dEta (same as 2D, but with 3D coordinates)
+ for (size_t i = 0; i < _nodes.size(); ++i)
+ jac.block(0, 0, 2, 3) += dsf.col(i) * _nodes[i]->getCoords().transpose();
+ // 2) compute thrid rows: dZeta = dXi X dEta
+ jac.row(2) = jac.row(0).cross(jac.row(1));
+ // 3) inverse Jacobian
+ ijac = jac.inverse();
+
+ // Get solution at element nodes
+ Eigen::VectorXd ue(_nodes.size());
+ for (size_t i = 0; i < _nodes.size(); ++i)
+ ue(i) = u[_nodes[i]->getId() - 1];
+
+ // Compute gradient by removing the third column of Jacobian
+ return ijac.block(0, 0, 3, 2) * dsf * ue;
}
diff --git a/sdpm/src/sdpmQuad4.h b/sdpm/src/sdpmQuad4.h
index 3f920db..4d9a948 100644
--- a/sdpm/src/sdpmQuad4.h
+++ b/sdpm/src/sdpmQuad4.h
@@ -29,20 +29,17 @@ namespace sdpm
*/
class SDPM_API Quad4 : public Element
{
-#ifndef SWIG
-protected:
- virtual void computeVolume() override;
- virtual void computeNormal() override;
-#endif
+ sdpmVector3d _tgt0; ///< unit tangent vector (node 0 to node 2)
+ sdpmVector3d _tgt1; ///< unit tangent vector (node 1 to node 3)
public:
- Quad4(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes) {}
- virtual ~Quad4() {}
+ Quad4(size_t n, Tag *tag, std::vector<Node *> &nodes);
+ ~Quad4() {}
virtual ElType getType() const override { return ElType::QUAD4; }
#ifndef SWIG
- virtual void init() override;
virtual void update() override;
+ virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const override;
#endif
};
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
new file mode 100644
index 0000000..ca8e0bd
--- /dev/null
+++ b/sdpm/src/sdpmSolver.cpp
@@ -0,0 +1,336 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpmSolver.h"
+#include "sdpmBuilder.h"
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include "sdpmMesh.h"
+#include "sdpmGmshExport.h"
+#include "sdpmResults.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+
+using namespace sdpm;
+
+Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
+{
+ // Say hi
+ std::cout << std::endl;
+ std::cout << "*******************************************************************************" << std::endl;
+ std::cout << "** Hi! My name is SDPM v1.0-2301 **" << std::endl;
+ std::cout << "** ULiege 2023-2023 **" << std::endl;
+ std::cout << "*******************************************************************************" << std::endl;
+
+ // Check the problem
+ pbl.check();
+
+ // Count panel, set up global->local id map and resize AIC matrices
+ _nP = 0;
+ int i = 0;
+ for (auto body : pbl.getBodies())
+ {
+ _nP += body->getElements().size();
+ for (auto e : body->getElements())
+ _rows[e] = i++;
+ }
+ _A = sdpmMatrixXd::Zero(_nP, _nP);
+ _B = sdpmMatrixXd::Zero(_nP, _nP);
+
+ // Set up variables
+ size_t ne = _pbl.getElements().size();
+ _phi.resize(ne, 0.);
+ _tau.resize(ne, 0.);
+ _mu.resize(ne, 0.);
+ _u.resize(ne, sdpmVector3d::Zero());
+ _cp.resize(ne, 0.);
+
+ // Display configuration
+ std::cout << "--- Physical groups ---\n";
+ std::cout << "Number of panels: " << _nP << "\n";
+ for (auto b : _pbl.getBodies())
+ std::cout << *b << "\n";
+ std::cout << "--- Freestream conditions ---\n"
+ << std::setprecision(5)
+ << "Angle of attack: " << _pbl.getAoa() * 180 / 3.14159 << " deg\n"
+ << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n";
+ std::cout << "--- Reference data ---\n"
+ << "Reference area: " << _pbl.getRefArea() << "\n"
+ << "Reference length: " << _pbl.getRefLgth() << "\n"
+ << "Reference origin: (" << _pbl.getRefCtr().transpose() << ")"
+ << "\n"
+ << std::endl;
+}
+
+Solver::~Solver()
+{
+ std::cout << "~sdpm::Solver()\n";
+}
+
+/**
+ * @brief Update elements and AIC matrices
+ */
+void Solver::update()
+{
+ // Update elements
+ _pbl.updateElements();
+
+ // Compute AIC
+ std::cout << "Computing AIC matrices... " << std::flush;
+ std::vector<Body *> bodies = _pbl.getBodies();
+ _A.setZero();
+ _B.setZero();
+ for (auto jbody : bodies)
+ {
+ // body
+ for (auto ej : jbody->getElements())
+ {
+ sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ej); // global to local transformation matrix
+ sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ej); // coordinates of source panel in local axes
+ for (auto ibody : bodies)
+ {
+ for (auto ei : ibody->getElements())
+ {
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l);
+ _A(_rows.at(ei), _rows.at(ej)) = mutau[0];
+ _B(_rows.at(ei), _rows.at(ej)) = mutau[1];
+ }
+ }
+ }
+ // wake
+ Wake const &wake = jbody->getWake();
+ std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+ for (auto ew : wake.getElements())
+ {
+ sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ew); // global to local transformation matrix
+ sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ew); // coordinates of source panel in local axes
+ for (auto ibody : bodies)
+ {
+ for (auto ei : ibody->getElements())
+ {
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l);
+ _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
+ _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+ }
+ }
+ }
+ }
+ // symmetry
+ if (_pbl.getSym())
+ {
+ for (auto jbody : bodies)
+ {
+ // body
+ for (auto ej : jbody->getElements())
+ {
+ sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ej); // global to local transformation matrix
+ sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ej); // coordinates of source panel in local axes
+ for (auto ibody : bodies)
+ {
+ for (auto ei : ibody->getElements())
+ {
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l);
+ _A(_rows.at(ei), _rows.at(ej)) -= mutau[0];
+ _B(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+ }
+ }
+ }
+ // wake
+ Wake const &wake = jbody->getWake();
+ std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+ for (auto ew : wake.getElements())
+ {
+ sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ew); // global to local transformation matrix
+ sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ew); // coordinates of source panel in local axes
+ for (auto ibody : bodies)
+ {
+ for (auto ei : ibody->getElements())
+ {
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l);
+ _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
+ _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+ }
+ }
+ }
+ }
+ }
+ std::cout << "done." << std::endl;
+}
+
+/**
+ * @brief Solve the potential equation
+ */
+void Solver::run()
+{
+ std::cout << "Running steady SDPM solver... " << std::endl;
+ std::cout << "computing sources... " << std::flush;
+ sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
+ sdpmVectorXd etau(_nP);
+ for (auto body : _pbl.getBodies())
+ for (auto e : body->getElements())
+ etau(_rows[e]) = ui.dot(e->getNormal());
+ std::cout << "done." << std::endl;
+
+ std::cout << "solving for doublets... " << std::flush;
+ sdpmVectorXd emu(_nP);
+ emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A).solve(_B * etau);
+ std::cout << "done." << std::endl;
+
+ std::cout << "computing flow and loads on bodies... " << std::flush;
+ computeFlow(etau, emu);
+ computeLoad();
+ std::cout << "done." << std::endl
+ << std::endl;
+}
+
+/**
+ * @brief Write the results
+ */
+void Solver::save(GmshExport const &writer, std::string const &suffix)
+{
+ // Write files
+ std::cout << "Saving files... " << std::endl;
+ // set up results
+ Results results;
+ results.scalarsAtElems["phi"] = &_phi;
+ results.scalarsAtElems["mu"] = &_mu;
+ results.scalarsAtElems["tau"] = &_tau;
+ results.vectorsAtElems["u"] = &_u;
+ results.scalarsAtElems["cp"] = &_cp;
+ // save results
+ writer.save(results, suffix);
+ // save aerodynamic loads breakdown
+ std::cout << "writing file: aeroloads_breakdown.dat... " << std::flush;
+ std::ofstream outfile;
+ outfile.open("aeroloads_breakdown.dat");
+ for (auto b : _pbl.getBodies())
+ {
+ outfile << "# Body - " << b->getName() << " (" << b->getElements().size() << " elements)" << std::endl;
+ outfile << std::fixed << std::setprecision(5)
+ << "Cl = " << b->cl << "\n"
+ << "Cd = " << b->cd << "\n"
+ << "Cs = " << b->cs << "\n"
+ << "Cm = " << b->cm << std::endl
+ << std::endl;
+ }
+ outfile.close();
+ std::cout << "done." << std::endl;
+}
+
+/**
+ * @brief Compute flow solution on bodies
+ */
+void Solver::computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
+{
+ // Store source and doublet strength and compute perturbation potential on elements
+ sdpmVectorXd ephi = -_A * emu - _B * etau;
+ std::vector<Body *> bodies = _pbl.getBodies();
+ for (auto body : bodies)
+ {
+ for (auto e : body->getElements())
+ {
+ size_t id = e->getId() - 1;
+ int i = _rows[e];
+ _tau[id] = etau(i);
+ _mu[id] = emu(i);
+ _phi[id] = ephi(i);
+ }
+ }
+
+ // Transfer doublet at nodes
+ std::vector<sdpmDouble> muNodes(_pbl.getNodes().size(), 0);
+ for (auto b : bodies)
+ {
+ // associate the value to its element
+ std::map<Element *, sdpmDouble> mapMu;
+ for (auto e : b->getElements())
+ mapMu[e] = _mu[e->getId() - 1];
+ // average at nodes
+ std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+ for (auto n : b->getNodes())
+ for (auto e : neMap.at(n))
+ muNodes[n->getId() - 1] += mapMu.at(e) / neMap.at(n).size();
+ }
+
+ // Compute surface flow functionals
+ sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
+ for (auto body : bodies)
+ {
+ for (auto e : body->getElements())
+ {
+ sdpmVector3d u = ui - _tau[e->getId() - 1] * e->getNormal() - e->computeGradient(muNodes);
+ _u[e->getId() - 1] = u; // velocity
+ _cp[e->getId() - 1] = 1 - u.dot(u); // pressure coefficient
+ }
+ }
+}
+
+/**
+ * @brief Compute aerodynamic load coefficients on bodies
+ */
+void Solver::computeLoad()
+{
+ sdpmDouble sRef = _pbl.getRefArea(); // reference surface area
+ sdpmDouble cRef = _pbl.getRefLgth(); // reference length
+ sdpmVector3d xRef = _pbl.getRefCtr(); // reference center
+ sdpmVector3d dragDir = _pbl.getDragDir(); // drag direction
+ sdpmVector3d sideDir = _pbl.getSideDir(); // sideforce direction
+ sdpmVector3d liftDir = _pbl.getLiftDir(); // lift direction
+ _cl = 0;
+ _cd = 0;
+ _cs = 0;
+ _cm = 0;
+ for (auto b : _pbl.getBodies())
+ {
+ // Reset
+ std::fill(b->nLoads.begin(), b->nLoads.end(), sdpmVector3d::Zero());
+
+ // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element and average at node
+ std::vector<Node *> nods = b->getNodes();
+ std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+ for (size_t i = 0; i < nods.size(); ++i)
+ for (auto e : neMap.at(nods[i]))
+ b->nLoads[i] += -_cp[e->getId() - 1] * e->getVolume() * e->getNormal() / e->getNodes().size();
+
+ // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
+ // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
+ b->cm = 0;
+ sdpmVector3d cf(0, 0, 0);
+ for (size_t i = 0; i < nods.size(); ++i)
+ {
+ sdpmVector3d l = nods[i]->getCoords() - xRef; // lever arm
+ sdpmVector3d cfi = b->nLoads[i]; // normalized force
+ cf += cfi;
+ b->cm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+ }
+ // rotate to flow direction
+ b->cd = cf.dot(dragDir) / sRef;
+ b->cs = cf.dot(sideDir) / sRef;
+ b->cl = cf.dot(liftDir) / sRef;
+ // compute total
+ _cl += b->cl;
+ _cd += b->cd;
+ _cs += b->cs;
+ _cm += b->cm;
+ }
+}
+
+void Solver::write(std::ostream &out) const
+{
+ out << "sdpm::Solver with\n"
+ << _pbl;
+}
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
new file mode 100644
index 0000000..16e0ebc
--- /dev/null
+++ b/sdpm/src/sdpmSolver.h
@@ -0,0 +1,76 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMSOLVER_H
+#define SDPMSOLVER_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include <vector>
+#include <map>
+
+namespace sdpm
+{
+
+/**
+ * @brief Source-doublet panel solver
+ * @authors Adrien Crovato
+ */
+class SDPM_API Solver : public sdpmObject
+{
+ // Problem
+ Problem &_pbl; ///< problem definition
+ // AICs
+ size_t _nP; ///< number of body panels
+ std::map<Element *, int> _rows; ///< map linking an element to a matrix cell
+ sdpmMatrixXd _A; ///< body-body doublet matrix
+ sdpmMatrixXd _B; ///< body-body source matrix
+ // Solution
+ std::vector<sdpmDouble> _phi; ///< perturbation potential
+ std::vector<sdpmDouble> _tau; ///< source strength
+ std::vector<sdpmDouble> _mu; ///< doublet strength
+ std::vector<sdpmVector3d> _u; ///< velocity
+ std::vector<sdpmDouble> _cp; ///< pressure coefficient
+ // Aerodynamic loads
+ sdpmDouble _cl; ///< lift coefficient
+ sdpmDouble _cd; ///< drag coefficient
+ sdpmDouble _cs; ///< sideforce coefficient
+ sdpmDouble _cm; ///< pitch moment coefficient (positive nose-up)
+
+public:
+ Solver(Problem &pbl);
+ ~Solver();
+
+ sdpmDouble getLiftCoeff() const { return _cl; }
+ sdpmDouble getDragCoeff() const { return _cd; }
+ sdpmDouble getSideforceCoeff() const { return _cs; }
+ sdpmDouble getMomentCoeff() const { return _cm; }
+
+ void update();
+ void run();
+ void save(GmshExport const &writer, std::string const &suffix = "");
+
+#ifndef SWIG
+ virtual void write(std::ostream &out) const override;
+#endif
+
+private:
+ void computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
+ void computeLoad();
+};
+
+} // namespace sdpm
+#endif // SDPMSOLVER_H
diff --git a/sdpm/src/sdpmTimers.h b/sdpm/src/sdpmTimers.h
index 6eae7a7..7786ec5 100644
--- a/sdpm/src/sdpmTimers.h
+++ b/sdpm/src/sdpmTimers.h
@@ -55,14 +55,14 @@ public:
* @brief Collection of timers
* @author Adrien Crovato
*/
-class SDPM_API Timers : public sdpm::sdpmObject
+class SDPM_API Timers : public sdpmObject
{
private:
std::map<std::string, Timer> _tms; ///< set of timers
public:
Timers() : sdpm::sdpmObject() {}
- virtual ~Timers() {}
+ ~Timers() {}
#ifndef SWIG
Timer &operator[](std::string const &name);
diff --git a/sdpm/src/sdpmUnitTest.cpp b/sdpm/src/sdpmUnitTest.cpp
index b239064..6f00cef 100644
--- a/sdpm/src/sdpmUnitTest.cpp
+++ b/sdpm/src/sdpmUnitTest.cpp
@@ -26,9 +26,13 @@ using namespace sdpm;
*/
void UnitTest::gmshIO(Mesh &mesh, GmshExport &writer)
{
+ // Get nodes and elements
+ std::vector<Node *> nods = mesh.getNodes();
+ std::vector<Element *> elems = mesh.getElements();
+
// Variables to be written
- std::vector<sdpmDouble> x(mesh.nodes.size(), 0.), type(mesh.elements.size(), 0.); // x-position of node, element type
- std::vector<sdpmVector3d> v(mesh.nodes.size(), sdpmVector3d::Zero()), cg(mesh.elements.size(), sdpmVector3d::Zero()); // x,y,z-position of node and element cg
+ std::vector<sdpmDouble> x(nods.size(), 0.), type(elems.size(), 0.); // x-position of node, element type
+ std::vector<sdpmVector3d> v(nods.size(), sdpmVector3d::Zero()), cg(elems.size(), sdpmVector3d::Zero()); // x,y,z-position of node and element cg
// Temporary results container
Results res;
@@ -38,16 +42,15 @@ void UnitTest::gmshIO(Mesh &mesh, GmshExport &writer)
res.vectorsAtElems["cg"] = &cg;
// Fill the variables and write them to disk
- for (size_t i = 0; i < mesh.nodes.size(); ++i)
+ for (size_t i = 0; i < nods.size(); ++i)
{
- x[i] = mesh.nodes[i]->getCoord()(0);
- v[i] = mesh.nodes[i]->getCoord();
+ x[i] = nods[i]->getCoords()(0);
+ v[i] = nods[i]->getCoords();
}
- for (size_t i = 0; i < mesh.elements.size(); ++i)
+ for (size_t i = 0; i < elems.size(); ++i)
{
- type[i] = static_cast<sdpmDouble>(mesh.elements[i]->getType());
- mesh.elements[i]->init();
- cg[i] = mesh.elements[i]->getCg();
+ type[i] = static_cast<sdpmDouble>(elems[i]->getType());
+ cg[i] = elems[i]->getCg();
}
- writer.save(mesh.name + "_new", res);
+ writer.save(res, "_new");
}
diff --git a/sdpm/src/sdpmUnitTest.h b/sdpm/src/sdpmUnitTest.h
index d8e24e9..77d5fe0 100644
--- a/sdpm/src/sdpmUnitTest.h
+++ b/sdpm/src/sdpmUnitTest.h
@@ -30,7 +30,7 @@ class SDPM_API UnitTest
{
public:
UnitTest() {}
- virtual ~UnitTest() {}
+ ~UnitTest() {}
void gmshIO(Mesh &mesh, GmshExport &writer);
};
diff --git a/sdpm/src/sdpmWake.cpp b/sdpm/src/sdpmWake.cpp
new file mode 100644
index 0000000..27a78aa
--- /dev/null
+++ b/sdpm/src/sdpmWake.cpp
@@ -0,0 +1,156 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#include "sdpmWake.h"
+#include "sdpmEdge.h"
+#include "sdpmNode.h"
+#include "sdpmElement.h"
+#include <vector>
+#include <unordered_set>
+using namespace sdpm;
+
+#include <iostream>
+
+/**
+ * @brief Create the wake in the case where the elements on the pressure side still have their original (suction side) trailing edge nodes
+ */
+Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag) : Group(msh, name)
+{
+ // Create TE edges
+ createEdges();
+
+ // Find wing elements having an edge on the TE
+ for (auto e : wingTag->getElements())
+ {
+ std::vector<Node *> nods = e->getNodes();
+ size_t idx = 0;
+ for (size_t i = 0; i < nods.size(); ++i)
+ {
+ // build wing edge and check if wing edge is a TE edge
+ std::vector<Node *> eN(2);
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = nods[(idx + j) % nods.size()];
+ Edge ed(eN);
+ auto it = _teEdges.find(&ed);
+ if (it != _teEdges.end())
+ {
+ // check normal direction
+ if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
+ (*it)->el1 = e; // suction side
+ else
+ (*it)->el2 = e; // pressure side
+ break;
+ }
+ idx++;
+ }
+ }
+
+ // Map wake to wing elements
+ createMap(wingTag);
+}
+
+/**
+ * @brief Create the wake in the case where the elements on the pressure side already have their correct (pressure side) trailing edge nodes
+ */
+Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap) : Group(msh, name)
+{
+ // Create TE edges
+ createEdges();
+
+ // Find wing elements having an edge on the TE
+ for (auto e : wingTag->getElements())
+ {
+ std::vector<Node *> nods = e->getNodes();
+ size_t idx = 0;
+ for (size_t i = 0; i < nods.size(); ++i, ++idx)
+ {
+ // build wing edge and check if wing edge is a TE edge
+ std::vector<Node *> eN(2);
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = nods[(idx + j) % nods.size()];
+ Edge ed(eN);
+ auto it = _teEdges.find(&ed);
+ // try with actual nodes
+ if (it != _teEdges.end())
+ {
+ (*it)->el1 = e; // suction side
+ break;
+ }
+ // otherwise try again, but by replacing the nodes by their corresponding upper nodes
+ else
+ {
+ // build wing edge and check if wing edge is a TE edge
+ std::vector<Node *> eN(2);
+ try
+ {
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = iTeMap.at(nods[(idx + j) % nods.size()]);
+ }
+ catch (const std::out_of_range &)
+ {
+ continue;
+ }
+ Edge ed(eN);
+ auto it = _teEdges.find(&ed);
+ if (it != _teEdges.end())
+ {
+ (*it)->el2 = e; // pressure side
+ break;
+ }
+ }
+ }
+ }
+
+ // Map wake to wing elements
+ createMap(wingTag);
+}
+
+/**
+ * @brief Create set of TE edges
+ */
+void Wake::createEdges()
+{
+ for (auto e : _tag->getElements())
+ {
+ std::vector<Node *> nods = e->getNodes();
+ Edge *ed = new Edge({nods[0], nods[3]});
+ ed->el0 = e;
+ _teEdges.insert(ed);
+ }
+}
+
+/**
+ * @brief Create a map between the wake and wing elements
+ */
+void Wake::createMap(Tag const *wingTag)
+{
+ for (auto ed : _teEdges)
+ {
+ if (ed->el1 == NULL || ed->el2 == NULL)
+ {
+ std::stringstream err;
+ err << "sdpm::Wake: at least one of the two trailing edge wing elements from tag " << *wingTag << "could not be associated to wake element " << *(ed->el0) << "!\n";
+ throw std::runtime_error(err.str());
+ }
+ _wwMap[ed->el0] = std::pair<Element *, Element *>(ed->el1, ed->el2);
+ delete ed;
+ }
+}
+
+void Wake::write(std::ostream &out) const
+{
+ out << "sdpm::Wake " << *_tag;
+}
diff --git a/sdpm/src/sdpmWake.h b/sdpm/src/sdpmWake.h
new file mode 100644
index 0000000..338892b
--- /dev/null
+++ b/sdpm/src/sdpmWake.h
@@ -0,0 +1,60 @@
+/*
+ * Copyright 2023 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.
+ */
+
+#ifndef SDPMWAKE_H
+#define SDPMWAKE_H
+
+#include "sdpm.h"
+#include "sdpmGroup.h"
+#include "sdpmTag.h"
+#include "sdpmEdge.h"
+#include <unordered_set>
+#include <map>
+
+namespace sdpm
+{
+
+/**
+ * @brief Hold the wake attached to a lifting body
+ * @authors Adrien Crovato
+ */
+class SDPM_API Wake : public Group
+{
+ std::unordered_set<Edge *, EquEdge, EquEdge> _teEdges; ///< trailing edge edges
+ std::map<Element *, std::pair<Element *, Element *>> _wwMap; ///< wing-wake map
+
+public:
+ Wake(Mesh &msh, std::string const &name, Tag *const &wingTag);
+ Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap);
+ ~Wake() {}
+
+private:
+ void createEdges();
+ void createMap(Tag const *tag);
+
+public:
+#ifndef SWIG
+ std::map<Element *, std::pair<Element *, Element *>> const &getMap() const
+ {
+ return _wwMap;
+ }
+ virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace sdpm
+
+#endif // SDPMWAKE_H
diff --git a/sdpm/src/sdpm_extras.h b/sdpm/src/sdpm_extras.h
index 55bbe5a..981800a 100644
--- a/sdpm/src/sdpm_extras.h
+++ b/sdpm/src/sdpm_extras.h
@@ -64,4 +64,4 @@ SDPM_API void tomatlab(std::string const &fname, Eigen::SparseMatrix<double> con
}
} // namespace sdpm
-#endif //SDPM_EXTRAS_H
+#endif // SDPM_EXTRAS_H
diff --git a/sdpm/tests/gmshio.py b/sdpm/tests/gmshio.py
index 817edc2..556f905 100644
--- a/sdpm/tests/gmshio.py
+++ b/sdpm/tests/gmshio.py
@@ -1,13 +1,13 @@
# -*- coding: utf-8 -*-
# Copyright 2023 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.
@@ -33,7 +33,7 @@ def main():
# write mesh
tms['msh_out'].start()
wrt = sdpm.GmshExport(msh)
- wrt.save(msh.name+'_new')
+ wrt.save('_new')
tst = sdpm.UnitTest()
tst.gmshIO(msh, wrt)
tms['msh_out'].stop()
@@ -41,13 +41,13 @@ def main():
print('Mesh info:')
print('-', msh)
print('Nodes info:')
- for n in msh.nodes:
+ for n in msh.getNodes():
print('-', n)
print('Elements info:')
- for e in msh.elements:
+ for e in msh.getElements():
print('-', e)
print('Tags info:')
- for n, t in msh.tags.items():
+ for n, t in msh.getTags().items():
print('-', n, '->', t)
print('')
# print timers
@@ -55,11 +55,11 @@ def main():
# test
tests = CTests()
- tests.add(CTest('Number of nodes', len(msh.nodes), 12, 1e-12))
- tests.add(CTest('Number of elements', len(msh.elements), 10, 1e-12))
- tests.add(CTest('Number of tags', len(msh.tags), 2, 1e-12))
+ tests.add(CTest('Number of nodes', len(msh.getNodes()), 12, 1e-12))
+ tests.add(CTest('Number of elements', len(msh.getElements()), 10, 1e-12))
+ tests.add(CTest('Number of tags', len(msh.getTags()), 2, 1e-12))
tests.run()
-
+
# eof
print('')
diff --git a/sdpm/tests/m6.py b/sdpm/tests/m6.py
new file mode 100644
index 0000000..eff0489
--- /dev/null
+++ b/sdpm/tests/m6.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+
+# 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.
+
+# Onera M6 wing
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import math
+
+def main():
+ # geometry parameters
+ pars = {'wNc': 100, 'wNs': 10, 'bC': 0.25, 'pS': 1.0}
+ # flow parameters
+ aoa = 3.06 * math.pi / 180
+
+ # start timer
+ tms = sdpm.Timers()
+ tms['total'].start()
+ # mesh
+ tms['msh'].start()
+ msh = MeshLoader(build_fpath('models/oneraM6.geo', __file__, 1)).run(pars)
+ tms['msh'].stop()
+ # problem
+ tms['pbl'].start()
+ pbl = sdpm.Problem(msh, aoa, 0, 0, 0.7528, 0.64, 0, 0, 0, True)
+ wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+ pbl.addBody(wing)
+ tms['pbl'].stop()
+ ## solver
+ tms['sol'].start()
+ sol = sdpm.Solver(pbl)
+ sol.update()
+ sol.run()
+ sol.save(sdpm.GmshExport(msh))
+ tms['sol'].stop()
+ # end timer
+ tms['total'].stop()
+ print(tms)
+
+ # test
+ tests = CTests()
+ tests.add(CTest('CL', sol.getLiftCoeff(), 0.200, 5e-2)) # gitlab.uliege.be/am-dept/dartflo: 0.200
+ tests.add(CTest('CD', sol.getDragCoeff(), 0.003, 5e-2)) # dartflo: 0.004
+ tests.add(CTest('CS', sol.getSideforceCoeff(), 0.008, 5e-2)) # dartflo: 0.011
+ tests.add(CTest('CM', sol.getMomentCoeff(), -0.147, 5e-2)) # dartflo: -0.148
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == '__main__':
+ main()
diff --git a/sdpm/tests/n12.py b/sdpm/tests/n12.py
new file mode 100644
index 0000000..d6fc9c4
--- /dev/null
+++ b/sdpm/tests/n12.py
@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2023 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.
+
+# Naca 0012 wing
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import math
+
+def main():
+ # geometry parameters
+ sym = 0
+ span = 5 # true span if sym=0, half span if sym=1
+ pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10}
+ # flow parameters
+ aoa = 3 * math.pi / 180
+
+ # start timer
+ tms = sdpm.Timers()
+ tms['total'].start()
+ # mesh
+ tms['msh'].start()
+ msh = MeshLoader(build_fpath('models/naca0012.geo', __file__, 1)).run(pars)
+ tms['msh'].stop()
+ # problem
+ tms['pbl'].start()
+ pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
+ wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+ pbl.addBody(wing)
+ tms['pbl'].stop()
+ ## solver
+ tms['sol'].start()
+ sol = sdpm.Solver(pbl)
+ sol.update()
+ sol.run()
+ sol.save(sdpm.GmshExport(msh))
+ tms['sol'].stop()
+ # end timer
+ tms['total'].stop()
+ print(tms)
+
+ # test
+ tests = CTests()
+ tests.add(CTest('CL', sol.getLiftCoeff(), 0.229, 5e-2)) # github.com/acrovato/aero: 0.226
+ tests.add(CTest('CD', sol.getDragCoeff(), 0.003, 5e-2)) # aero: 0.003
+ tests.add(CTest('CS', sol.getSideforceCoeff(), 0.000, 5e-2))
+ tests.add(CTest('CM', sol.getMomentCoeff(), -0.056, 5e-2))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == '__main__':
+ main()
diff --git a/sdpm/tests/n12tail.py b/sdpm/tests/n12tail.py
new file mode 100644
index 0000000..aee9152
--- /dev/null
+++ b/sdpm/tests/n12tail.py
@@ -0,0 +1,77 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2023 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.
+
+# Naca 0012 wing and tail
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import math
+
+def main():
+ # geometry parameters
+ sym = 0
+ span = 5 # wing (true span if sym=0, half span if sym=1)
+ tspan = 3 # tail (true span if sym=0, half span if sym=1)
+ pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10, 'tSpn': tspan, 'tNc': 100, 'tNs': 6}
+ # flow parameters
+ aoa = 3 * math.pi / 180
+
+ # start timer
+ tms = sdpm.Timers()
+ tms['total'].start()
+ # mesh
+ tms['msh'].start()
+ msh = MeshLoader(build_fpath('models/naca0012.geo', __file__, 1)).run(pars)
+ tms['msh'].stop()
+ # problem
+ tms['pbl'].start()
+ pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
+ wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+ tail = sdpm.Body(msh, 'tail', 'tTe', 15)
+ pbl.addBody(wing)
+ pbl.addBody(tail)
+ tms['pbl'].stop()
+ ## solver
+ tms['sol'].start()
+ sol = sdpm.Solver(pbl)
+ sol.update()
+ sol.run()
+ sol.save(sdpm.GmshExport(msh))
+ tms['sol'].stop()
+ # end timer
+ tms['total'].stop()
+ print(tms)
+
+ # test
+ tests = CTests()
+ tests.add(CTest('Wing CL', wing.cl, 0.230, 5e-2))
+ tests.add(CTest('Wing CD', wing.cd, 0.003, 5e-2))
+ tests.add(CTest('Wing CS', wing.cs, 0.000, 5e-2))
+ tests.add(CTest('Wing CM', wing.cm, -0.056, 5e-2))
+ tests.add(CTest('Tail CL', tail.cl, 0.045, 5e-2))
+ tests.add(CTest('Tail CD', tail.cd, 0.001, 5e-2))
+ tests.add(CTest('Tail CS', tail.cs, 0.000, 5e-2))
+ tests.add(CTest('Tail CM', tail.cm, -0.229, 5e-2))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == '__main__':
+ main()
--
GitLab