diff --git a/sdpm/CMakeLists.txt b/sdpm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..dfa6b125378eee97b52e57c4df1aa2277ae1fada
--- /dev/null
+++ b/sdpm/CMakeLists.txt
@@ -0,0 +1,28 @@
+# 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.
+
+# Add source dir
+ADD_SUBDIRECTORY( src )
+ADD_SUBDIRECTORY( _src )
+
+# Add test dir
+MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
+
+# Add to install
+INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
+ ${CMAKE_CURRENT_LIST_DIR}/utils.py
+ DESTINATION sdpm)
+INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
+ ${CMAKE_CURRENT_LIST_DIR}/viscous
+ DESTINATION sdpm)
diff --git a/sdpm/__init__.py b/sdpm/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..b2401ae7a4f7757f92a1d10ba2ec8efc6022e36f
--- /dev/null
+++ b/sdpm/__init__.py
@@ -0,0 +1,21 @@
+# -*- coding: utf-8 -*-
+# dart MODULE initialization file
+
+# 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.
+
+
+import fwk
+import tbox
+from sdpmw import *
diff --git a/sdpm/_src/CMakeLists.txt b/sdpm/_src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..38adbaead554a94fa178f13d9426e07a46ea102f
--- /dev/null
+++ b/sdpm/_src/CMakeLists.txt
@@ -0,0 +1,35 @@
+# 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.
+
+# CMake input file of the SWIG wrapper around "spmw.so"
+
+INCLUDE(${SWIG_USE_FILE})
+
+FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
+FILE(GLOB ISRCS *.i)
+
+SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS "-interface" "_sdpmw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(sdpmw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET sdpmw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(sdpmw)
+
+TARGET_INCLUDE_DIRECTORIES(sdpmw PRIVATE ${PROJECT_SOURCE_DIR}/sdpm/_src
+ ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
+ ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
+ ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(sdpmw PRIVATE sdpm tbox fwk ${PYTHON_LIBRARIES})
+
+INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/sdpmw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS sdpmw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/sdpm/_src/sdpmw.h b/sdpm/_src/sdpmw.h
new file mode 100644
index 0000000000000000000000000000000000000000..ee069ab530e48f8f14def0dbb60a31af1f6f755d
--- /dev/null
+++ b/sdpm/_src/sdpmw.h
@@ -0,0 +1,18 @@
+/*
+ * 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.
+ */
+
+#include "pPanels.h"
+#include "pSDSolver.h"
\ No newline at end of file
diff --git a/sdpm/_src/sdpmw.i b/sdpm/_src/sdpmw.i
new file mode 100644
index 0000000000000000000000000000000000000000..e039abc6192354cff08b699bcc1e4b2b1b858f14
--- /dev/null
+++ b/sdpm/_src/sdpmw.i
@@ -0,0 +1,47 @@
+/*
+ * 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.
+ */
+
+// SWIG input file of the 'dart' module
+
+%feature("autodoc","1");
+
+%module(docstring=
+"'sdpmw' module: Viscous inviscid interaction for dart 2021-2022
+(c) ULiege - A&M",
+directors="0",
+threads="1"
+) sdpmw
+%{
+
+#include "sdpm.h"
+
+#include "fwkw.h"
+#include "tboxw.h"
+#include "sdpmw.h"
+
+%}
+
+
+%include "fwkw.swg"
+
+// ----------- MODULES UTILISES ------------
+%import "tboxw.i"
+
+// ----------- FLOW CLASSES ----------------
+%include "sdpm.h"
+
+%include "pPanels.h"
+%include "pSDSolver.h"
\ No newline at end of file
diff --git a/sdpm/src/CMakeLists.txt b/sdpm/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..aa9c304e587744d5df3bbef761c8a07129b5ec32
--- /dev/null
+++ b/sdpm/src/CMakeLists.txt
@@ -0,0 +1,27 @@
+# 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.
+
+# CMake input file of dart.so
+
+FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
+
+ADD_LIBRARY(sdpm SHARED ${SRCS})
+MACRO_DebugPostfix(sdpm)
+TARGET_INCLUDE_DIRECTORIES(sdpm PUBLIC ${PROJECT_SOURCE_DIR}/sdpm/src)
+
+TARGET_LINK_LIBRARIES(sdpm tbox)
+
+INSTALL(TARGETS sdpm DESTINATION ${CMAKE_INSTALL_PREFIX})
+
+SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
diff --git a/sdpm/src/pPanels.cpp b/sdpm/src/pPanels.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9305723db3f61994fdbd445c33126f7c1d6f8b7
--- /dev/null
+++ b/sdpm/src/pPanels.cpp
@@ -0,0 +1,657 @@
+#include "pPanels.h"
+#include<iostream>
+#include <Eigen/Dense>
+#include <cmath>
+#include "math.h"
+
+using namespace sdpm;
+
+Panels::Panels(size_t _nPanChord)
+{
+ nPanChord = _nPanChord;
+ nPanSpan = 2;
+ beta = 1;
+}
+
+void Panels::CreateWing()
+{
+ Eigen::RowVectorXd xp(2*nPanChord+1); // Non-linear chordwise distribution
+ Eigen::RowVectorXd xc(2*nPanChord+1); // Control point chordwise distribution
+
+ Eigen::RowVectorXd dumm = Eigen::VectorXd::LinSpaced(2*nPanChord+1, 0, M_PI);
+ for (auto i=0; i<xp.size(); ++i)
+ {
+ xp(i) = 1-std::sin(dumm(i));
+ }
+
+ for (auto i=0; i<xp.size(); ++i)
+ {
+ xc(i) = (xp(i+1) + xp(i))/2;
+ }
+
+ /* Non-linear spanwise distribution */
+
+ Eigen::RowVectorXd yp(nPanSpan+1);
+ Eigen::RowVectorXd yc(nPanSpan);
+ Eigen::RowVectorXd theta = Eigen::RowVectorXd::LinSpaced(nPanSpan+1, M_PI, 0);
+ for (auto i=0; i<theta.size(); ++i)
+ {
+ yp(i) = std::cos(theta(i));
+ }
+
+ for (auto i=0; i<yc.size(); ++i)
+ {
+ yc(i) = (yp(i+1) + yp(i)) / 2;
+ }
+
+ /* Spanwise chord variation */
+ Eigen::RowVectorXd cvar;
+ cvar = (lambda-1)*c0*yp.cwiseAbs();
+ cvar += c0 * Eigen::MatrixXd::Ones(1, cvar.size());
+ Eigen::RowVectorXd epsilonvar;
+ epsilonvar=epsilon*yp.cwiseAbs();
+
+ /* Spanwise quarter chord distance due to twist */
+ Eigen::RowVectorXd xsweep;
+ xsweep = std::tan(LambdaC4)*yp.cwiseAbs()*b/2;
+ /* Spanwise vertical distance due to dihedral */
+ Eigen::RowVectorXd zdih;
+ zdih = std::tan(dihedral)*yp.cwiseAbs()*b/2;
+
+
+ /* Calculate airfoil shape at each panel gridpoint */
+
+ Eigen::MatrixXd zp = Eigen::MatrixXd::Zero(1,2*nPanChord+1);
+
+ double t = .12; // Airfoil maximum thickness
+
+ /* Upper side */
+
+ for (auto i=0; i<zp.size()/2; ++i)
+ {
+ zp(i) = -(t/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
+ }
+
+ /* Lower side */
+
+ for (auto i=zp.size()/2; i<zp.size(); ++i)
+ {
+ zp(i) = (t/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
+ }
+
+ // iko=find(xp <= maxpos);
+ // zp(iko)=zp(iko)+maxcamb/maxpos^2*(2*maxpos*xp(iko)-xp(iko).^2);
+ // iko=find(xp > maxpos);
+ // zp(iko)=zp(iko)+maxcamb/(1-maxpos)^2*((1-2*maxpos)+2*maxpos*xp(iko)-xp(iko).^2);
+
+ /* Assemble bound panel corner point matrices */
+ /* Eigen::Array * Eigen::Array is a mult element wise while Eigen::Matrix * Eigen::Matrix is a matrix product -> cwiseProduct */
+ Xp = xp.transpose().replicate(1,nPanSpan+1).cwiseProduct(cvar.replicate(2*nPanChord+1,1)) - cvar.replicate(2*nPanChord+1,1)*ax;
+ Yp = yp.replicate(2*nPanChord+1,1) * b/2;
+ Zp = zp.transpose().replicate(1,nPanSpan+1).cwiseProduct(cvar.replicate(2*nPanChord+1,1));
+
+ Surface = (c0+c0*lambda)/2 * b;
+
+
+
+ // % Apply twist
+ // for i=1:n+1
+ // Ry=(cos(epsilonvar(i)) sin(epsilonvar(i));-sin(epsilonvar(i)) cos(epsilonvar(i)));
+ // dummy=Ry*(Xp(:,i)';Zp(:,i)');
+ // Xp(:,i)=dummy(1,:)';
+ // Zp(:,i)=dummy(2,:)';
+ // end
+
+ // % Apply dihedral
+ // Zp=Zp+repmat(zdih,2*m+1,1);
+ // % Apply sweep
+ // Xp=Xp+repmat(xsweep,2*m+1,1);
+
+ Xc.setZero(2*nPanChord, nPanSpan);
+ Yc.setZero(2*nPanChord, nPanSpan);
+ Zc.setZero(2*nPanChord, nPanSpan);
+ nx.setZero(2*nPanChord, nPanSpan);
+ ny.setZero(2*nPanChord, nPanSpan);
+ nz.setZero(2*nPanChord, nPanSpan);
+ Xc0.setZero(2*nPanChord, nPanSpan);
+ Yc0.setZero(2*nPanChord, nPanSpan);
+ Zc0.setZero(2*nPanChord, nPanSpan);
+ nx0.setZero(2*nPanChord, nPanSpan);
+ ny0.setZero(2*nPanChord, nPanSpan);
+ nz0.setZero(2*nPanChord, nPanSpan);
+ tauxx.setZero(2*nPanChord, nPanSpan);
+ tauxy.setZero(2*nPanChord, nPanSpan);
+ tauxz.setZero(2*nPanChord, nPanSpan);
+ tauyx.setZero(2*nPanChord, nPanSpan);
+ tauyy.setZero(2*nPanChord, nPanSpan);
+ tauyz.setZero(2*nPanChord, nPanSpan);
+ s.setZero(2*nPanChord, nPanSpan);
+ s0.setZero(2*nPanChord, nPanSpan);
+ cpln.setZero(2*nPanChord, nPanSpan);
+
+ for (size_t i=0; i<2*nPanChord; ++i)
+ {
+ for (size_t j=0; j<nPanSpan; ++j)
+ {
+ /* Calculate control points */
+ Xc0(i,j) = (Xp(i,j) + Xp(i,j+1) + Xp(i+1,j+1) + Xp(i+1,j))/4;
+ Yc0(i,j) = (Yp(i,j) + Yp(i,j+1) + Yp(i+1,j+1) + Yp(i+1,j))/4;
+ Zc0(i,j) = (Zp(i,j) + Zp(i,j+1) + Zp(i+1,j+1) + Zp(i+1,j))/4;
+
+ /* Calculate the normal vectors */
+ /*(N,s0(i,j),cpln(i,j))=normalarea((Xp(i,j) Xp(i,j+1) Xp(i+1,j+1) Xp(i+1,j)), ...
+ (Yp(i,j) Yp(i,j+1) Yp(i+1,j+1) Yp(i+1,j)), ...
+ (Zp(i,j) Zp(i,j+1) Zp(i+1,j+1) Zp(i+1,j))); */
+
+ /* Normal vector (cross product) */
+ Eigen::RowVector3d normal1 = Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i,j+1)-Xp(i+1,j), Yp(i,j+1)-Yp(i+1,j), Zp(i,j+1)-Zp(i+1,j)));
+ Eigen::RowVectorXd mag1 = (normal1*normal1.transpose());
+ mag1(0) = std::sqrt(mag1(0));
+ normal1 /= mag1(0); // mag1 is of size 1
+
+ /* Calculate a second normal from points 1,2,3 only */
+ Eigen::RowVector3d normal2= Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i,j+1)-Xp(i,j), Yp(i,j+1)-Yp(i,j), Zp(i,j+1)-Zp(i,j)));
+ Eigen::RowVectorXd mag2=(normal2*normal2.transpose());
+ mag2(0) = std::sqrt(mag2(0));
+ normal2=normal2/mag2(0);
+
+ /* Test if the two normals are colinear */
+ cpln(i,j) = normal1.cross(normal2)(0);
+ Eigen::RowVectorXd dum1 = cpln*cpln.transpose();
+ cpln(i,j)=sqrt(dum1(0));
+
+ /* Calculate a third normal from points 1,3,4 only */
+ Eigen::RowVector3d normal3 = Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i+1,j)-Xp(i,j), Yp(i+1,j)-Yp(i,j), Zp(i+1,j)-Zp(i,j)));
+ Eigen::RowVectorXd mag3=(normal3*normal3.transpose());
+ mag3(0) = std::sqrt(mag3(0));
+ /* Panel area */
+ s0(i,j) = mag2(0)/2+mag3(0)/2;
+
+ nx0(i,j) = normal1(0);
+ ny0(i,j) = normal1(1);
+ nz0(i,j) = normal1(2);
+
+ }
+ }
+
+ /* Scale for compressibility */
+ /*Xp0=Xp;
+ Yp0=Yp;
+ Zp0=Zp;
+ Xp=Xp0/beta; // (beta = sqrt(1-M*M))
+ Yp=Yp0;
+ Zp=Zp0;*/
+
+ for (size_t i=0 ; i<2*nPanChord; ++i)
+ {
+ for (size_t j=0; j<nPanSpan; ++j)
+ {
+ /* Calculate control points */
+ Xc(i,j)=(Xp(i,j) + Xp(i,j+1) + Xp(i+1,j+1) + Xp(i+1,j))/4;
+ Yc(i,j)=(Yp(i,j) + Yp(i,j+1) + Yp(i+1,j+1) + Yp(i+1,j))/4;
+ Zc(i,j)=(Zp(i,j) + Zp(i,j+1) + Zp(i+1,j+1) + Zp(i+1,j))/4;
+ /* Calculate the normal vectors */
+ /* (N,s(i,j),cpln(i,j))=normalarea((Xp(i,j) Xp(i,j+1) Xp(i+1,j+1) Xp(i+1,j)), ...
+ (Yp(i,j) Yp(i,j+1) Yp(i+1,j+1) Yp(i+1,j)), ...
+ (Zp(i,j) Zp(i,j+1) Zp(i+1,j+1) Zp(i+1,j))); */
+
+ /* Normal vector (cross product) */
+ Eigen::RowVector3d normal1 = Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i,j+1)-Xp(i+1,j), Yp(i,j+1)-Yp(i+1,j), Zp(i,j+1)-Zp(i+1,j)));
+ Eigen::RowVectorXd mag1 = (normal1*normal1.transpose());
+ mag1(0) = std::sqrt(mag1(0));
+ normal1 /= mag1(0); // mag1 is of size 1
+
+ /* Calculate a second normal from points 1,2,3 only */
+ Eigen::RowVector3d normal2= Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i,j+1)-Xp(i,j), Yp(i,j+1)-Yp(i,j), Zp(i,j+1)-Zp(i,j)));
+ Eigen::RowVectorXd mag2=(normal2*normal2.transpose());
+ mag2(0) = std::sqrt(mag2(0));
+ normal2=normal2/mag2(0);
+
+ /* Test if the two normals are colinear */
+ cpln(i,j) = normal1.cross(normal2)(0);
+ Eigen::RowVectorXd dum1 = cpln*cpln.transpose();
+ cpln(i,j)=sqrt(dum1(0));
+
+ /* Calculate a third normal from points 1,3,4 only */
+ Eigen::RowVector3d normal3 = Eigen::RowVector3d(Xp(i+1,j+1)-Xp(i,j), Yp(i+1,j+1)-Yp(i,j), Zp(i+1,j+1)-Zp(i,j)).cross(Eigen::RowVector3d(Xp(i+1,j)-Xp(i,j), Yp(i+1,j)-Yp(i,j), Zp(i+1,j)-Zp(i,j)));
+ Eigen::RowVectorXd mag3=(normal3*normal3.transpose());
+ mag3(0) = std::sqrt(mag3(0));
+ /* Panel area */
+ s(i,j) = mag2(0)/2+mag3(0)/2;
+
+ nx(i,j) = normal1(0);
+ ny(i,j) = normal1(1);
+ nz(i,j) = normal1(2);
+
+ /* Calculate tangential vector in x direction */
+ double a=1/sqrt(1+(normal1(0)/normal1(2))*(normal1(0)/normal1(2)));
+ double c=-normal1(0)/normal1(2)*a;
+ tauxx(i,j)=a;
+ tauxy(i,j)=0;
+ tauxz(i,j)=c;
+
+ /* Calculate tangential vector in y direction */
+ Eigen::RowVector3d dummy = normal1.cross(Eigen::RowVector3d(a, 0, c));
+ tauyx(i,j)=dummy(0);
+ tauyy(i,j)=dummy(1);
+ tauyz(i,j)=dummy(2);
+ }
+ }
+
+}
+
+void Panels::WakeControlPoints(unsigned long it){
+ for (size_t i=0; i<it; ++i){
+ for (size_t j=0; j<nPanSpan; ++j)
+ {
+
+ /* Calculate control points */
+ Xcw(i,j) = (Xw(i,j) + Xw(i,j+1) + Xw(i+1,j+1) + Xw(i+1,j))/4;
+ Ycw(i,j) = (Yw(i,j) + Yw(i,j+1) + Yw(i+1,j+1) + Yw(i+1,j))/4;
+ Zcw(i,j) = (Zw(i,j) + Zw(i,j+1) + Zw(i+1,j+1) + Zw(i+1,j))/4;
+ /* Calculate the normal vectors */
+ /* N=normalarea((Xw(i,j) Xw(i,j+1) Xw(i+1,j+1) Xw(i+1,j)), ...
+ (Yw(i,j) Yw(i,j+1) Yw(i+1,j+1) Yw(i+1,j)), ...
+ (Zw(i,j) Zw(i,j+1) Zw(i+1,j+1) Zw(i+1,j))); */
+
+ /* Normal vector (cross product) */
+ Eigen::RowVector3d normal1 = Eigen::RowVector3d(Xw(i+1,j+1)-Xw(i,j), Yw(i+1,j+1)-Yw(i,j), Zw(i+1,j+1)-Zw(i,j)).cross(Eigen::RowVector3d(Xw(i,j+1)-Xw(i+1,j), Yw(i,j+1)-Yw(i+1,j), Zw(i,j+1)-Zw(i+1,j)));
+ Eigen::RowVectorXd mag1 = (normal1*normal1.transpose());
+ mag1(0) = std::sqrt(mag1(0));
+ normal1 /= mag1(0); // mag1 is of size 1
+
+ nxw(i,j) = normal1(0);
+ nyw(i,j) = normal1(1);
+ nzw(i,j) = normal1(2);
+
+ /* Calculate tangential vector in x direction */
+ double a= 1/sqrt(1+(normal1(0)/normal1(2))*(normal1(0)/normal1(2)));
+ double c=-normal1(0)/normal1(2)*a;
+ tauxxw(i,j)=a;
+ tauxyw(i,j)=0;
+ tauxzw(i,j)=c;
+
+ /* Calculate tangential vector in y direction */
+ Eigen::Vector3d dummy = normal1.cross(Eigen::Vector3d(a, 0, c));
+ tauyxw(i,j)=dummy(0);
+ tauyyw(i,j)=dummy(1);
+ tauyzw(i,j)=dummy(2);
+ }
+ }
+}
+
+
+void Panels::ComputeInfluenceCoeffs_wing()
+{
+ /* Sources */
+ AN.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Au.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Av.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Aw.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ /* Doublets */
+ BN.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Bu.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Bv.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ Bw.setZero(2*nPanChord*nPanSpan,2*nPanChord*nPanSpan);
+ /* Call mex function to calculate influence coefficients of wing on wing */
+
+ Eigen::RowVectorXd x, y, z;
+ x.setZero(5);
+ y.setZero(5);
+ z.setZero(5);
+
+ Eigen::RowVector3d normal;
+ Eigen::RowVector3d taux;
+ Eigen::RowVector3d tauy;
+
+ Eigen::RowVector2d Phisd;
+ Eigen::RowVectorXd uvw;
+ uvw.setZero(6);
+
+ double x0, y0, z0, xcp, ycp, zcp;
+
+ size_t ic,jc,i,j;
+
+ size_t m = Xc.rows();
+ size_t n = Xc.cols();
+ //size_t mp1=m+1;
+ //size_t np1=n+1;
+ //size_t mn=m*n;
+
+ /* Create Influence Matrix */
+
+ for (ic=0; ic < m; ic++) {
+ for (jc=0; jc < n; jc++) {
+ for (i=0; i < m; i++) {
+ for (j=0; j < n; j++) {
+ x0=Xc(ic,jc);
+ y0=Yc(ic,jc);
+ z0=Zc(ic,jc);
+ xcp=Xc(i,j);
+ ycp=Yc(i,j);
+ zcp=Zc(i,j);
+ x(0)=Xp(i,j);
+ x(1)=Xp(i,(j+1));
+ x(2)=Xp(i+1,(j+1));
+ x(3)=Xp(i+1,j);
+ x(4)=Xp(i,j);
+ y(0)=Yp(i,j);
+ y(1)=Yp(i,(j+1));
+ y(2)=Yp(i+1,(j+1));
+ y(3)=Yp(i+1,j);
+ y(4)=Yp(i,j);
+ z(0)=Zp(i,j);
+ z(1)=Zp(i,(j+1));
+ z(2)=Zp(i+1,(j+1));
+ z(3)=Zp(i+1,j);
+ z(4)=Zp(i,j);
+ normal(0)=nx(i,j);
+ normal(1)=ny(i,j);
+ normal(2)=nz(i,j);
+ taux(0)=tauxx(i,j);
+ taux(1)=tauxy(i,j);
+ taux(2)=tauxz(i,j);
+ tauy(0)=tauyx(i,j);
+ tauy(1)=tauyy(i,j);
+ tauy(2)=tauyz(i,j);
+
+ sdpanel_wing(Phisd,uvw,x0,y0,z0,x,y,z,xcp,ycp,zcp,normal,taux,tauy);
+ AN(ic*n+jc,i*n+j)=Phisd(0);
+ Bu(ic*n+jc,i*n+j)=-uvw(3);
+ Bv(ic*n+jc,i*n+j)=-uvw(4);
+ Bw(ic*n+jc,i*n+j)=-uvw(5);
+ if ((ic*n+jc) == (i*n+j)) {
+ Au(ic*n+jc,i*n+j)=normal(0)/2.0;
+ Av(ic*n+jc,i*n+j)=normal(1)/2.0;
+ Aw(ic*n+jc,i*n+j)=normal(2)/2.0;
+ BN(ic*n+jc,i*n+j)=-0.5;
+ }else{
+ Au(ic*n+jc,i*n+j)=uvw(0);
+ Av(ic*n+jc,i*n+j)=uvw(1);
+ Aw(ic*n+jc,i*n+j)=uvw(2);
+ BN(ic*n+jc,i*n+j)=Phisd(1);
+ }
+
+ }
+ }
+ }
+ }
+}
+
+void Panels::sdpanel_wing(Eigen::RowVector2d &Phisd, Eigen::RowVectorXd &uvw, double x0, double y0, double z0, Eigen::RowVectorXd x, Eigen::RowVectorXd y, Eigen::RowVectorXd z, double xcp, double ycp, double zcp,
+ Eigen::RowVector3d normal, Eigen::RowVector3d taux, Eigen::RowVector3d tauy)
+{
+
+ Eigen::RowVectorXd xl, yl, zl, r, e, h, mij, dij, rij, wij;
+ xl.setZero(5);
+ yl.setZero(5);
+ zl.setZero(5);
+ r.setZero(5);
+ e.setZero(5);
+ h.setZero(5);
+
+ mij.setZero(4);
+ dij.setZero(4);
+ rij.setZero(4);
+ wij.setZero(4);
+
+ double x0l,y0l,z0l;
+ double pi,Phi1,Phi2,logterm,logdenom,atanterm,fraction,fracdenom,xx,yy,us,vs,ws,ud,vd,wd;
+ int i;
+
+ pi=3.14159265358979;
+
+ for (i=0; i < 5; i++) {
+ xl(i)=(x(i)-xcp)*taux(0)+(y(i)-ycp)*taux(1)+(z(i)-zcp)*taux(2);
+ yl(i)=(x(i)-xcp)*tauy(0)+(y(i)-ycp)*tauy(1)+(z(i)-zcp)*tauy(2);
+ zl(i)=(x(i)-xcp)*normal(0)+(y(i)-ycp)*normal(1)+(z(i)-zcp)*normal(2);
+ }
+ x0l=(x0-xcp)*taux(0)+(y0-ycp)*taux(1)+(z0-zcp)*taux(2);
+ y0l=(x0-xcp)*tauy(0)+(y0-ycp)*tauy(1)+(z0-zcp)*tauy(2);
+ z0l=(x0-xcp)*normal(0)+(y0-ycp)*normal(1)+(z0-zcp)*normal(2);
+
+ for (i=0; i < 4; i++) {
+ r(i)=sqrt((x0l-xl(i))*(x0l-xl(i))+(y0l-yl(i))*(y0l-yl(i))+z0l*z0l);
+ rij(i)=(x0l-xl(i))*(x0l-xl(i+1))+(y0l-yl(i))*(y0l-yl(i+1))+z0l*z0l;
+ wij(i)=(x0l-xl(i+1))*(y0l-yl(i))-(x0l-xl(i))*(y0l-yl(i+1));
+ e(i)=(x0l-xl(i))*(x0l-xl(i))+z0l*z0l;
+ h(i)=(x0l-xl(i))*(y0l-yl(i));
+ mij(i)=(yl(i+1)-yl(i))/(xl(i+1)-xl(i));
+ dij(i)=sqrt((xl(i+1)-xl(i))*(xl(i+1)-xl(i))+(yl(i+1)-yl(i))*(yl(i+1)-yl(i)));
+ }
+ r(4)=r(0);
+ e(4)=e(0);
+ h(4)=h(0);
+
+ Phi1=0.;
+ Phi2=0.;
+ us=0.;
+ vs=0.;
+ ud=0.;
+ vd=0.;
+ wd=0.;
+ for (i=0; i < 4; i++) {
+ logterm=0.0;
+ logdenom=(r(i)+r(i+1)-dij(i));
+ if (fabs(logdenom) > 1e-10){
+ logterm=log((r(i)+r(i+1)+dij(i))/logdenom);
+ }
+ us=us-(yl(i+1)-yl(i))*logterm/dij(i);
+ vs=vs-(xl(i)-xl(i+1))*logterm/dij(i);
+ Phi1+=((x0l-xl(i))*(yl(i+1)-yl(i))-(y0l-yl(i))*(xl(i+1)-xl(i)))/dij(i)*logterm;
+ atanterm=0.0;
+ if ((fabs(xl(i+1)-xl(i)) > 1e-10) && (fabs(z0l) > 1e-10)){ //mij(i) is not infinite
+ xx=(mij(i)*e(i)-h(i))/(z0l*r(i));
+ yy=(mij(i)*e(i+1)-h(i+1))/(z0l*r(i+1));
+ atanterm=atan2(xx-yy,1+xx*yy);
+ }
+ fraction=0.0;
+ fracdenom=r(i)*r(i+1)*(r(i)*r(i+1)+rij(i));
+ if (fabs(fracdenom) > 1e-10){
+ fraction=(r(i)+r(i+1))/fracdenom;
+ }
+ Phi2+=atanterm;
+ ud+=z0l*(yl(i)-yl(i+1))*fraction;
+ vd+=z0l*(xl(i+1)-xl(i))*fraction;
+ wd+=wij(i)*fraction;
+ }
+ // Source potential
+ Phisd(0)=-1/4.0/pi*(Phi1+fabs(z0l)*Phi2);
+// Phisd(0)=-1/4.0/pi*(Phi1-z0l*Phi2);
+ // Source velocities
+ us=us/4.0/pi;
+ vs=vs/4.0/pi;
+ ws=Phi2/4.0/pi;
+ // Doublet potential
+ Phisd(1)=ws;
+ // Doublet velocities
+ ud=ud/4.0/pi;
+ vd=vd/4.0/pi;
+ wd=wd/4.0/pi;
+
+ // Transform velocities to global coordinates
+ uvw(0)=taux(0)*us+tauy(0)*vs+normal(0)*ws;
+ uvw(1)=taux(1)*us+tauy(1)*vs+normal(1)*ws;
+ uvw(2)=taux(2)*us+tauy(2)*vs+normal(2)*ws;
+ uvw(3)=taux(0)*ud+tauy(0)*vd+normal(0)*wd;
+ uvw(4)=taux(1)*ud+tauy(1)*vd+normal(1)*wd;
+ uvw(5)=taux(2)*ud+tauy(2)*vd+normal(2)*wd;
+}
+
+void Panels::ComputeInfluenceCoeffs_wake(unsigned long it)
+{
+
+ CN.setZero(2*nPanChord*nPanSpan, nPanSpan*it);
+ Cu.setZero(2*nPanChord*nPanSpan, nPanSpan*it);
+ Cv.setZero(2*nPanChord*nPanSpan, nPanSpan*it);
+ Cw.setZero(2*nPanChord*nPanSpan, nPanSpan*it);
+
+ Eigen::RowVectorXd x, y, z;
+ x.setZero(5);
+ y.setZero(5);
+ z.setZero(5);
+
+ Eigen::RowVector3d normal;
+ Eigen::RowVector3d taux;
+ Eigen::RowVector3d tauy;
+
+ Eigen::RowVectorXd Phisd;
+ Phisd.setZero(1);
+ Eigen::RowVectorXd uvw;
+ uvw.setZero(3);
+
+ double x0, y0, z0, xcp, ycp, zcp;
+
+ size_t ic,jc,i,j;
+
+ size_t m = Xc.rows();
+ size_t n = Xc.cols();
+ //size_t mp1=m+1;
+ //size_t np1=n+1;
+ //size_t mn=m*n;
+
+ size_t mwp1 = it+1;
+ size_t mw = mwp1-1;
+
+ /* Create Influence Matrix */
+
+ for (ic=0; ic < m; ic++) {
+ for (jc=0; jc < n; jc++) {
+ for (i=0; i < mw; i++) {
+ for (j=0; j < n; j++) {
+ x0=Xc(ic,jc);
+ y0=Yc(ic,jc);
+ z0=Zc(ic,jc);
+ xcp=Xcw(i,j);
+ ycp=Ycw(i,j);
+ zcp=Zcw(i,j);
+ x(0)=Xw(i,j);
+ x(1)=Xw(i,(j+1));
+ x(2)=Xw(i+1,(j+1));
+ x(3)=Xw(i+1,j);
+ x(4)=Xw(i,j);
+ y(0)=Yw(i,j);
+ y(1)=Yw(i,(j+1));
+ y(2)=Yw(i+1,(j+1));
+ y(3)=Yw(i+1,j);
+ y(4)=Yw(i,j);
+ z(0)=Zw(i,j);
+ z(1)=Zw(i,(j+1));
+ z(2)=Zw(i+1,(j+1));
+ z(3)=Zw(i+1,j);
+ z(4)=Zw(i,j);
+ normal(0)=nxw(i,j);
+ normal(1)=nyw(i,j);
+ normal(2)=nzw(i,j);
+ taux(0)=tauxxw(i,j);
+ taux(1)=tauxyw(i,j);
+ taux(2)=tauxzw(i,j);
+ tauy(0)=tauyxw(i,j);
+ tauy(1)=tauyyw(i,j);
+ tauy(2)=tauyzw(i,j);
+
+ sdpanel_wake(Phisd,uvw,x0,y0,z0,x,y,z,xcp,ycp,zcp,normal,taux,tauy);
+ CN(ic*n+jc, i*n+j)=Phisd(0);
+ Cu(ic*n+jc, i*n+j)=-uvw(0);
+ Cv(ic*n+jc, i*n+j)=-uvw(1);
+ Cw(ic*n+jc, i*n+j)=-uvw(2);
+ }
+ }
+ }
+ }
+}
+
+void Panels::sdpanel_wake(Eigen::RowVectorXd &Phisd, Eigen::RowVectorXd &uvw,double x0, double y0, double z0, Eigen::RowVectorXd x, Eigen::RowVectorXd y, Eigen::RowVectorXd z, double xcp, double ycp, double zcp,
+ Eigen::RowVector3d normal, Eigen::RowVector3d taux, Eigen::RowVector3d tauy)
+{
+
+ Eigen::RowVectorXd xl, yl, zl, r, e, h, mij, dij, rij, wij;
+ xl.setZero(5);
+ yl.setZero(5);
+ zl.setZero(5);
+ r.setZero(5);
+ e.setZero(5);
+ h.setZero(5);
+
+ mij.setZero(4);
+ dij.setZero(4);
+ rij.setZero(4);
+ wij.setZero(4);
+
+ double x0l,y0l,z0l;
+ double pi,Phi2,atanterm,xx,yy,fraction,fracdenom,ud,vd,wd;
+ int i;
+
+ pi=3.14159265358979;
+
+ for (i=0; i < 5; i++) {
+ xl(i)=(x(i)-xcp)*taux(0)+(y(i)-ycp)*taux(1)+(z(i)-zcp)*taux(2);
+ yl(i)=(x(i)-xcp)*tauy(0)+(y(i)-ycp)*tauy(1)+(z(i)-zcp)*tauy(2);
+ zl(i)=(x(i)-xcp)*normal(0)+(y(i)-ycp)*normal(1)+(z(i)-zcp)*normal(2);
+ }
+ x0l=(x0-xcp)*taux(0)+(y0-ycp)*taux(1)+(z0-zcp)*taux(2);
+ y0l=(x0-xcp)*tauy(0)+(y0-ycp)*tauy(1)+(z0-zcp)*tauy(2);
+ z0l=(x0-xcp)*normal(0)+(y0-ycp)*normal(1)+(z0-zcp)*normal(2);
+
+ for (i=0; i < 4; i++) {
+ r(i)=sqrt((x0l-xl(i))*(x0l-xl(i))+(y0l-yl(i))*(y0l-yl(i))+z0l*z0l);
+ rij(i)=(x0l-xl(i))*(x0l-xl(i+1))+(y0l-yl(i))*(y0l-yl(i+1))+z0l*z0l;
+ wij(i)=(x0l-xl(i+1))*(y0l-yl(i))-(x0l-xl(i))*(y0l-yl(i+1));
+ e(i)=(x0l-xl(i))*(x0l-xl(i))+z0l*z0l;
+ h(i)=(x0l-xl(i))*(y0l-yl(i));
+ mij(i)=(yl(i+1)-yl(i))/(xl(i+1)-xl(i));
+ dij(i)=sqrt((xl(i+1)-xl(i))*(xl(i+1)-xl(i))+(yl(i+1)-yl(i))*(yl(i+1)-yl(i)));
+ }
+ r(4)=r(0);
+ e(4)=e(0);
+ h(4)=h(0);
+
+ Phi2=0;
+ ud=0.0;
+ vd=0.0;
+ wd=0.0;
+ for (i=0; i < 4; i++) {
+ atanterm=0;
+ if ((xl(i+1)-xl(i) != 0) && (z0l != 0)){ //mij(i) is not infinite
+ xx=(mij(i)*e(i)-h(i))/(z0l*r(i));
+ yy=(mij(i)*e(i+1)-h(i+1))/(z0l*r(i+1));
+ atanterm=atan2(xx-yy,1+xx*yy);
+ }
+ fraction=0.0;
+ fracdenom=r(i)*r(i+1)*(r(i)*r(i+1)+rij(i));
+ if (fabs(fracdenom) > 1e-10){
+ fraction=(r(i)+r(i+1))/fracdenom;
+ }
+ Phi2+=atanterm;
+ ud+=z0l*(yl(i)-yl(i+1))*fraction;
+ vd+=z0l*(xl(i+1)-xl(i))*fraction;
+ wd+=wij(i)*fraction;
+ }
+ // Doublet potential
+ Phisd(0)=Phi2/4.0/pi;
+ // Doublet velocities
+ ud=ud/4.0/pi;
+ vd=vd/4.0/pi;
+ wd=wd/4.0/pi;
+
+ // Transform velocities to global coordinates
+ uvw(0)=taux(0)*ud+tauy(0)*vd+normal(0)*wd;
+ uvw(1)=taux(1)*ud+tauy(1)*vd+normal(1)*wd;
+ uvw(2)=taux(2)*ud+tauy(2)*vd+normal(2)*wd;
+}
+
+
+
+void Panels::SetKutta()
+{
+ /* Set up Kutta condition */
+ Eigen::MatrixXd Kutta_dummy;
+ Kutta_dummy.setZero(nPanSpan,2*nPanChord*nPanSpan);
+ for (size_t i=0; i<nPanSpan; ++i)
+ {
+ Kutta_dummy(i,i)=1;
+ Kutta_dummy(i,(2*nPanChord-1)*nPanSpan+i)=-1;
+ }
+ Kutta.setZero(Kutta_dummy.rows(), Kutta_dummy.cols()+(Eigen::MatrixXd::Identity(nPanSpan, nPanSpan).cols()));
+ Kutta << Kutta_dummy, Eigen::MatrixXd::Identity(nPanSpan, nPanSpan);
+}
\ No newline at end of file
diff --git a/sdpm/src/pPanels.h b/sdpm/src/pPanels.h
new file mode 100644
index 0000000000000000000000000000000000000000..736330cf28afda60ab2a5784812da3e01cd6ffb8
--- /dev/null
+++ b/sdpm/src/pPanels.h
@@ -0,0 +1,135 @@
+#ifndef PPANELS_H
+#define PPANELS_H
+
+#include "sdpm.h"
+#include <vector>
+#include <Eigen/Dense>
+
+namespace sdpm
+{
+
+class DART_API Panels
+{
+
+private:
+
+public:
+ double c0 = 1; /* Root chord (m) */
+ double lambda = 1; /* Taper ratio (-) */
+ double epsilon = 0; /* Twist angle (rad) */
+ double LambdaC4 = 0; /* Sweep angle at 1/4 chord (rad) */
+ double b = 1; /* Span (m) */
+ double dihedral = 0; /* Dihedral angle (rad) */
+ double ax = 0; /* Axis around which twist and taper are defined */
+
+ double beta;
+ double Surface;
+
+ size_t nPanChord, nPanSpan;
+
+ /* Bound panel corner point matrices */
+ Eigen::MatrixXd Xp; // Size (2m+1, n)
+ Eigen::MatrixXd Yp;
+ Eigen::MatrixXd Zp;
+
+ /* Compressibility scaling */
+ Eigen::MatrixXd Xp0; // =Xp;
+ Eigen::MatrixXd Yp0; // =Yp;
+ Eigen::MatrixXd Zp0; // =Zp;
+
+ /* Control points */
+ Eigen::MatrixXd Xc; // =zeros(2*m,n);
+ Eigen::MatrixXd Yc;
+ Eigen::MatrixXd Zc;
+ Eigen::MatrixXd nx;
+ Eigen::MatrixXd ny;
+ Eigen::MatrixXd nz;
+ Eigen::MatrixXd tauxx;
+ Eigen::MatrixXd tauxy;
+ Eigen::MatrixXd tauxz;
+ Eigen::MatrixXd tauyx;
+ Eigen::MatrixXd tauyy;
+ Eigen::MatrixXd tauyz;
+
+ Eigen::MatrixXd Xc0;
+ Eigen::MatrixXd Yc0;
+ Eigen::MatrixXd Zc0;
+
+ Eigen::MatrixXd s0;
+ Eigen::MatrixXd cpln;
+
+ /* Surface */
+ Eigen::MatrixXd s;
+
+ /* Wake */
+ Eigen::MatrixXd Xw;
+ Eigen::MatrixXd Yw;
+ Eigen::MatrixXd Zw;
+
+ /* Wake */
+ Eigen::MatrixXd MAw; // (mw, n)
+ Eigen::MatrixXd Xcw;
+ Eigen::MatrixXd Ycw;
+ Eigen::MatrixXd Zcw;
+ Eigen::MatrixXd nxw;
+ Eigen::MatrixXd nyw;
+ Eigen::MatrixXd nzw;
+ Eigen::MatrixXd tauxxw;
+ Eigen::MatrixXd tauxyw;
+ Eigen::MatrixXd tauxzw;
+ Eigen::MatrixXd tauyxw;
+ Eigen::MatrixXd tauyyw;
+ Eigen::MatrixXd tauyzw;
+ Eigen::MatrixXd cp; // (2mw, n)
+
+ /* Normal vectors of compressible geometry panels */
+ Eigen::MatrixXd nx0;
+ Eigen::MatrixXd ny0;
+ Eigen::MatrixXd nz0;
+
+ /* Influence coefficients */
+
+ /* Sources */
+ Eigen::MatrixXd AN;
+ Eigen::MatrixXd Au;
+ Eigen::MatrixXd Av;
+ Eigen::MatrixXd Aw;
+ /* Doublets */
+ Eigen::MatrixXd BN;
+ Eigen::MatrixXd Bu;
+ Eigen::MatrixXd Bv;
+ Eigen::MatrixXd Bw;
+
+ /* Wake Doublets only */
+ Eigen::MatrixXd CN;
+ Eigen::MatrixXd Cu;
+ Eigen::MatrixXd Cv;
+ Eigen::MatrixXd Cw;
+
+ Eigen::MatrixXd Kutta;
+
+ Panels(size_t nPanChord);
+
+ //void CreateWake(double U, double V, double W, double dt);
+ void ComputeInfluenceCoeffs_wing();
+ void ComputeInfluenceCoeffs_wake(unsigned long it);
+
+ void sdpanel_wing(Eigen::RowVector2d &Phisd, Eigen::RowVectorXd &uvw, double x0, double y0, double z0, Eigen::RowVectorXd x, Eigen::RowVectorXd y, Eigen::RowVectorXd z, double xcp, double ycp, double zcp,
+ Eigen::RowVector3d normal, Eigen::RowVector3d taux, Eigen::RowVector3d tauy);
+
+ void sdpanel_wake(Eigen::RowVectorXd &Phisd, Eigen::RowVectorXd &uvw,double x0, double y0, double z0, Eigen::RowVectorXd x, Eigen::RowVectorXd y, Eigen::RowVectorXd z, double xcp, double ycp, double zcp,
+ Eigen::RowVector3d normal, Eigen::RowVector3d taux, Eigen::RowVector3d tauy);
+
+ void WakeControlPoints(unsigned long it);
+
+
+ void SetKutta();
+
+
+ double Getc0(){return c0;};
+
+ void CreateWing();
+
+};
+} // namespace dart
+#endif //PPANELS_H
diff --git a/sdpm/src/pSDSolver.cpp b/sdpm/src/pSDSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..721fc449736eb65f33348697a18a8887193c6211
--- /dev/null
+++ b/sdpm/src/pSDSolver.cpp
@@ -0,0 +1,770 @@
+#include "pSDSolver.h"
+#include "pPanels.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+using namespace sdpm;
+
+SDSolver::SDSolver(Panels *_pan, double _alpha, double _airspeed, unsigned long _maxIt)
+{
+ pan = _pan;
+ fs_alpha = _alpha;
+ fs_airspeed = _airspeed;
+ maxIt = _maxIt;
+}
+
+void SDSolver::Run()
+{
+
+ int R = 1; /* Wake type */
+
+ Eigen::Vector3d fs_velocity;
+
+ fs_velocity(0) = fs_airspeed * std::cos(fs_alpha);
+ fs_velocity(1) = 0;
+ fs_velocity(2) = fs_airspeed * std::sin(fs_alpha);
+
+ double dt = pan->c0/pan->nPanChord/fs_velocity(0);
+
+ pan->Xw.setZero(maxIt+1, pan->nPanSpan+1);
+ pan->Yw.setZero(maxIt+1, pan->nPanSpan+1);
+ pan->Zw.setZero(maxIt+1, pan->nPanSpan+1);
+
+ Eigen::MatrixXd MAw = Eigen::MatrixXd::Zero(maxIt, pan->nPanSpan);
+ pan->Xcw.setZero(maxIt, pan->nPanSpan);
+ pan->Ycw.setZero(maxIt, pan->nPanSpan);
+ pan->Zcw.setZero(maxIt, pan->nPanSpan);
+ pan->nxw.setZero(maxIt, pan->nPanSpan);
+ pan->nyw.setZero(maxIt, pan->nPanSpan);
+ pan->nzw.setZero(maxIt, pan->nPanSpan);
+ pan->tauxxw.setZero(maxIt, pan->nPanSpan);
+ pan->tauxyw.setZero(maxIt, pan->nPanSpan);
+ pan->tauxzw.setZero(maxIt, pan->nPanSpan);
+ pan->tauyxw.setZero(maxIt, pan->nPanSpan);
+ pan->tauyyw.setZero(maxIt, pan->nPanSpan);
+ pan->tauyzw.setZero(maxIt, pan->nPanSpan);
+
+ pan->ComputeInfluenceCoeffs_wing();
+
+ pan->Xp0 = pan->Xp;
+ pan->Yp0 = pan->Yp;
+ pan->Zp0 = pan->Zp;
+
+ /* std::cout << " xp yp zp " << std::endl;
+ std::cout << pan->Xp << std::endl;
+ std::cout << pan->Yp << std::endl;
+ std::cout << pan->Zp << std::endl;
+
+ std::cout << " xc yc zc " << std::endl;
+ std::cout << pan->Xc << std::endl;
+ std::cout << pan->Yc << std::endl;
+ std::cout << pan->Zc << std::endl; */
+
+ /* Calculate derivative of doublet strength in the chordwise direction */
+ /* dX */
+ Eigen::MatrixXd dXm1 = (pan->Xc.row(1)-pan->Xc.row(0));
+ Eigen::MatrixXd dXm2 = (pan->Xc.block(2,0,2*pan->nPanChord,pan->Xc.cols()) - pan->Xc.block(0,0,2*pan->nPanChord-2,pan->Xc.cols()));
+ Eigen::MatrixXd dXm3 = (pan->Xc.row(2*pan->nPanChord-1)-pan->Xc.row(2*pan->nPanChord-2));
+
+ /* Vertical concatenation */
+ Eigen::MatrixXd dXm(dXm1.rows()+dXm2.rows()+dXm3.rows(), dXm1.cols());
+ dXm << dXm1,
+ dXm2,
+ dXm3;
+
+ /* dY */
+ Eigen::MatrixXd dYm1 = (pan->Yc.row(1)-pan->Yc.row(0));
+ Eigen::MatrixXd dYm2 = (pan->Yc.block(2,0,2*pan->nPanChord,pan->Yc.cols()) - pan->Yc.block(0,0,2*pan->nPanChord-2,pan->Yc.cols()));
+ Eigen::MatrixXd dYm3 = (pan->Yc.row(2*pan->nPanChord-1)-pan->Yc.row(2*pan->nPanChord-2));
+
+ /* Vertical concatenation */
+ Eigen::MatrixXd dYm(dYm1.rows()+dYm2.rows()+dYm3.rows(), dYm1.cols());
+ dYm << dYm1,
+ dYm2,
+ dYm3;
+
+ /* dZ */
+ Eigen::MatrixXd dZm1 = (pan->Zc.row(1)-pan->Zc.row(0));
+ Eigen::MatrixXd dZm2 = (pan->Zc.block(2,0,2*pan->nPanChord,pan->Zc.cols()) - pan->Zc.block(0,0,2*pan->nPanChord-2,pan->Zc.cols()));
+ Eigen::MatrixXd dZm3 = (pan->Zc.row(2*pan->nPanChord-1)-pan->Zc.row(2*pan->nPanChord-2));
+
+ /* Vertical concatenation */
+ Eigen::MatrixXd dZm(dZm1.rows()+dZm2.rows()+dZm3.rows(), dZm1.cols());
+ dZm << dZm1,
+ dZm2,
+ dZm3;
+
+ Eigen::MatrixXd sm = (dXm.array()*dXm.array() + dYm.array()*dYm.array() + dZm.array()).sqrt();
+
+ /* Chordwise tangent vector */
+ Eigen::MatrixXd tmx = dXm.array()/sm.array();
+ Eigen::MatrixXd tmy = dYm.array()/sm.array();
+ Eigen::MatrixXd tmz = dZm.array()/sm.array();
+
+ /* std::cout << " tm " << std::endl;
+ std::cout << tmx << std::endl;
+ std::cout << tmy << std::endl;
+ std::cout << tmz << std::endl; */
+
+
+ Eigen::MatrixXd dXn;
+ Eigen::MatrixXd dYn;
+ Eigen::MatrixXd dZn;
+
+ if (pan->nPanSpan > 2)
+ {
+ /* Calculate derivative of doublet strength in the spanwise direction */
+ /* dX */
+ Eigen::MatrixXd dXn1 = (pan->Xc.col(1)-pan->Xc.col(0));
+ Eigen::MatrixXd dXn2 = (pan->Xc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Xc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ Eigen::MatrixXd dXn3 = (pan->Xc.col(pan->nPanSpan-1)-pan->Xc.col(pan->nPanSpan-2));
+
+ /* Horizontal concatenation */
+ dXn.resize(dXn1.rows(), dXn1.cols()+dXn2.cols()+dXn3.cols());
+ dXn << dXn1, dXn2, dXn3;
+
+ /* dY */
+ Eigen::MatrixXd dYn1 = (pan->Yc.col(1)-pan->Yc.col(0));
+ Eigen::MatrixXd dYn2 = (pan->Yc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Yc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ Eigen::MatrixXd dYn3 = (pan->Yc.col(pan->nPanSpan-1)-pan->Yc.col(pan->nPanSpan-2));
+
+ /* Horizontal concatenation */
+ dYn.resize(dYn1.rows(), dYn1.cols()+dYn2.cols()+dYn3.cols());
+ dYn << dYn1, dYn2, dYn3;
+
+ /* dZ */
+ Eigen::MatrixXd dZn1 = (pan->Zc.col(1)-pan->Zc.col(0));
+ Eigen::MatrixXd dZn2 = (pan->Zc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Zc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ Eigen::MatrixXd dZn3 = (pan->Zc.col(pan->nPanSpan-1)-pan->Zc.col(pan->nPanSpan-2));
+
+ /* Horizontal concatenation */
+ dZn.resize(dZn1.rows(), dZn1.cols()+dZn2.cols()+dZn3.cols());
+ dZn << dZn1, dZn2, dZn3;
+
+ }
+ else if (pan->nPanSpan == 2){
+
+ /* Calculate derivative of doublet strength in the spanwise direction */
+ /* dX */
+ Eigen::MatrixXd dXn1 = (pan->Xc.col(1)-pan->Xc.col(0));
+
+ /* Horizontal concatenation */
+ dXn.resize(dXn1.rows(), dXn1.cols()+dXn1.cols());
+ dXn << dXn1, dXn1;
+
+ /* dY */
+ Eigen::MatrixXd dYn1 = (pan->Yc.col(1)-pan->Yc.col(0));
+
+ /* Horizontal concatenation */
+ dYn.resize(dYn1.rows(), dYn1.cols()+dYn1.cols());
+ dYn << dYn1, dYn1;
+
+ /* dZ */
+ Eigen::MatrixXd dZn1 = (pan->Zc.col(1)-pan->Zc.col(0));
+
+ /* Horizontal concatenation */
+ dZn.resize(dZn1.rows(), dZn1.cols()+dZn1.cols());
+ dZn << dZn1, dZn1;
+
+ }
+ else{
+ throw;
+ }
+
+ Eigen::MatrixXd sn = (dXn.array()*dXn.array() + dYn.array()*dYn.array() + dZn.array()).sqrt();
+
+ /* std::cout << " dxm dym dzm " << std::endl;
+ std::cout << dXm << std::endl;
+ std::cout << dYm << std::endl;
+ std::cout << dZm << std::endl;
+
+ std::cout << " dxn dyn dzn " << std::endl;
+ std::cout << dXn << std::endl;
+ std::cout << dYn << std::endl;
+ std::cout << dZn << std::endl; */
+
+
+ /* Calculate spanwise tangent vector */
+
+ Eigen::MatrixXd tnx = dXn.array()/sn.array();
+ Eigen::MatrixXd tny = dYn.array()/sn.array();
+ Eigen::MatrixXd tnz = dZn.array()/sn.array();
+
+ /* std::cout << " tnx tny tnz " << std::endl;
+ std::cout << tnx << std::endl;
+ std::cout << tny << std::endl;
+ std::cout << tnz << std::endl; */
+
+ /* Calculate source strengths */
+
+ Eigen::MatrixXd SA = -(fs_velocity[0]*pan->nx + fs_velocity[1]*pan->ny + fs_velocity[2]*pan->nz);
+ Eigen::MatrixXd dum = SA.transpose();
+ Eigen::MatrixXd sigma = dum.reshaped(2*pan->nPanSpan*pan->nPanChord,1);
+
+ /* std::cout << " SA sigma " << std::endl;
+ std::cout << SA << std::endl;
+ std::cout << sigma << std::endl; */
+
+ /* Solution at zero time */
+
+ Eigen::MatrixXd LHS = pan->BN;
+ Eigen::MatrixXd RHS = -pan->AN * sigma;
+ Eigen::VectorXd Sol = LHS.partialPivLu().solve(RHS);
+
+ Eigen::MatrixXd Phi0 = Sol.reshaped(pan->nPanSpan, 2*pan->nPanChord).transpose();
+
+ /* std::cout << "Phi0 " << Phi0 << std::endl; */
+
+ /* First wake row */
+ pan->Xw.row(0) = pan->Xp.row(2*pan->nPanChord);
+ pan->Yw.row(0) = pan->Yp.row(2*pan->nPanChord);
+ pan->Zw.row(0) = pan->Zp.row(2*pan->nPanChord);
+
+ CL.setZero(maxIt);
+ CD.setZero(maxIt);
+ CY.setZero(maxIt);
+
+ std::vector<Eigen::MatrixXd> Phi(maxIt);
+
+ unsigned long it = 1;
+
+ while(it < maxIt){
+
+ std::cout << it << std::endl;
+
+ /* Propagate wake */
+ Eigen::MatrixXd uwp = Eigen::MatrixXd::Zero(it-1, pan->nPanSpan + 1);
+ Eigen::MatrixXd vwp = Eigen::MatrixXd::Zero(it-1, pan->nPanSpan + 1);
+ Eigen::MatrixXd wwp = Eigen::MatrixXd::Zero(it-1, pan->nPanSpan + 1);
+
+ Eigen::MatrixXd dum0(1+uwp.rows(), uwp.cols());
+ dum0 << Eigen::MatrixXd::Zero(1,pan->nPanSpan+1),
+ uwp;
+
+ pan->Xw.block(1,0,it,pan->Xw.cols()) = pan->Xw.block(0,0,it-1,pan->Xw.cols()) + (Eigen::MatrixXd::Constant(dum0.rows(), dum0.cols(), fs_velocity(0)) + dum0)*dt;
+
+ Eigen::MatrixXd dum1(1+vwp.rows(), vwp.cols());
+ dum1 << Eigen::MatrixXd::Zero(1,pan->nPanSpan+1),
+ vwp;
+ pan->Yw.block(1,0,it,pan->Yw.cols()) = pan->Yw.block(0,0,it-1,pan->Yw.cols()) + (Eigen::MatrixXd::Constant(dum1.rows(), dum1.cols(), fs_velocity(1)) + dum1)*dt;
+
+ Eigen::MatrixXd dum2(1+uwp.rows(), uwp.cols());
+ dum2 << Eigen::MatrixXd::Zero(1,pan->nPanSpan+1),
+ wwp;
+ pan->Zw.block(1,0,it,pan->Zw.cols()) = pan->Zw.block(0,0,it-1,pan->Zw.cols()) + (Eigen::MatrixXd::Constant(dum2.rows(), dum2.cols(), fs_velocity(2)) + dum2)*dt;
+
+ if (it > 1){
+ MAw.block(1,0,it ,MAw.cols()) = MAw.block(0,0,it-1,MAw.cols());
+ }
+
+ pan->Xw.row(0) = pan->Xp.row(2*pan->nPanChord);
+ pan->Yw.row(0) = pan->Yp.row(2*pan->nPanChord);
+ pan->Zw.row(0) = pan->Zp.row(2*pan->nPanChord);
+
+ /* Calculate wake control points */
+
+ pan->WakeControlPoints(it);
+
+ /* Calculate wake influence coefficients */
+
+ pan->ComputeInfluenceCoeffs_wake(it);
+
+ /* Set up Kutta condition */
+
+ pan->SetKutta();
+
+ /* Equations */
+
+ Eigen::MatrixXd LHS;
+
+ /* std::cout << " Here " << std::endl;
+ std::cout << pan->CN.block(0,0,pan->CN.rows(), pan->nPanSpan) << std::endl; */
+
+ /* Horizontal concatenation */
+ Eigen::MatrixXd LHS_dum1(pan->BN.rows(), pan->BN.cols()+pan->CN.block(0,0,pan->CN.rows(), pan->nPanSpan).cols());
+ LHS_dum1 << pan->BN, pan->CN.block(0,0,pan->CN.rows(), pan->nPanSpan);
+
+ /* Vertical concatenation */
+ LHS.setZero(LHS_dum1.rows()+pan->Kutta.rows(), LHS_dum1.cols());
+ LHS << LHS_dum1,
+ pan->Kutta;
+
+ Eigen::MatrixXd RHS;
+
+ if (it==1){
+ Eigen::MatrixXd RHS_dum1;
+ RHS_dum1 = -pan->AN*sigma;
+
+ RHS.setZero(RHS_dum1.rows()+Eigen::MatrixXd::Zero(pan->nPanSpan, 1).rows(), RHS_dum1.cols());
+ RHS << RHS_dum1,
+ Eigen::MatrixXd::Zero(pan->nPanSpan, 1);
+ }
+ else{
+
+ Eigen::MatrixXd RHS_dum1;
+ RHS_dum1 = -pan->AN*sigma - pan->CN.block(0,pan->nPanSpan, pan->CN.rows(), (it-1)*pan->nPanSpan) * MAw.block(1,0,it-1,MAw.cols()).transpose().reshaped((it-1)*pan->nPanSpan,1);
+
+ RHS.setZero(RHS_dum1.rows()+Eigen::MatrixXd::Zero(pan->nPanSpan, 1).rows(), RHS_dum1.cols());
+ RHS << RHS_dum1,
+ Eigen::MatrixXd::Zero(pan->nPanSpan, 1);
+
+ }
+
+ /* Solve for strength */
+
+ /* std::cout << " LHS " << std::endl;
+ std::cout << LHS << std::endl;
+
+ std::cout << " RHS " << std::endl;
+ std::cout << RHS << std::endl; */
+
+ Eigen::VectorXd Mu2;
+ Mu2.setZero(RHS.size());
+ Mu2 = LHS.partialPivLu().solve(RHS);
+
+ /* std::cout << " Mu2 " << std::endl;
+ std::cout << Mu2 << std::endl; */
+
+ /*std::cout << "mu2" << std::endl;
+ std::cout << Mu2 << std::endl; */
+
+ Eigen::MatrixXd MA_dum1 = Mu2.reshaped(pan->nPanSpan, 2*pan->nPanChord+1).transpose();
+
+ /* std::cout << "madum" << std::endl;
+ std::cout << MA_dum1 << std::endl; */
+
+ MAw.row(0) = MA_dum1.row(2*pan->nPanChord);
+
+ /* std::cout << "maw" << std::endl;
+ std::cout << MAw << std::endl; */
+
+ Eigen::VectorXd Mu(2*pan->nPanChord*pan->nPanSpan);
+
+
+ for (size_t i=0; i<2*pan->nPanChord*pan->nPanSpan; ++i)
+ {
+ Mu(i) = Mu2(i);
+ }
+ /* std::cout << "mu" << std::endl;
+ std::cout << Mu << std::endl; */
+
+ Eigen::MatrixXd MA(2*pan->nPanChord, MA_dum1.cols());
+
+ for (size_t i=0; i<2*pan->nPanChord; ++i)
+ {
+ MA.row(i) = MA_dum1.row(i);
+ }
+
+ /* std::cout << "ma" << std::endl;
+ std::cout << MA << std::endl; */
+
+ Phi[it] = MA;
+
+ /* std::cout << " Phi[it] " << Phi[it] << std::endl; */
+
+ /* std::cout << " phi it " << std::endl;
+ std::cout << Phi[it] << std::endl; */
+
+ /* dMA_chord */
+ Eigen::MatrixXd dMA1_chord = (MA.row(1)-MA.row(0));
+ Eigen::MatrixXd dMA2_chord = (MA.block(2,0,2*pan->nPanChord-1,MA.cols()) - MA.block(0,0,2*pan->nPanChord-2,MA.cols()));
+ Eigen::MatrixXd dMA3_chord = (MA.row(2*pan->nPanChord-1)-MA.row(2*pan->nPanChord-2));
+
+ /* Vertical concatenation */
+ Eigen::MatrixXd dMA_chord(dMA1_chord.rows()+dMA2_chord.rows()+dMA3_chord.rows(), dMA1_chord.cols());
+ dMA_chord << dMA1_chord,
+ dMA2_chord,
+ dMA3_chord;
+
+ Eigen::MatrixXd um = dMA_chord.array()/sm.array();
+
+ /* dMA_span */
+
+ Eigen::MatrixXd dMA_span;
+
+ if (pan->nPanSpan > 2){
+ Eigen::MatrixXd dMA1_span = (MA.col(1)-MA.col(0));
+ Eigen::MatrixXd dMA2_span = (MA.block(0,2,2*pan->nPanChord,pan->nPanSpan) - MA.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ Eigen::MatrixXd dMA3_span = (MA.col(pan->nPanSpan-1)-MA.col(pan->nPanSpan-2));
+
+ /* Horizontal concatenation */
+ dMA_span.setZero(dMA1_span.rows(), dMA1_span.cols()+dMA2_span.cols()+dMA3_span.cols());
+ dMA_span << dMA1_span, dMA2_span, dMA3_span;
+ }
+ else if (pan->nPanSpan == 2){
+ Eigen::MatrixXd dMA1_span = (MA.col(1)-MA.col(0));
+
+ /* Horizontal concatenation */
+ dMA_span.setZero(dMA1_span.rows(), dMA1_span.cols()+dMA1_span.cols());
+ dMA_span << dMA1_span, dMA1_span;
+ }
+ else{
+ throw;
+ }
+
+ Eigen::MatrixXd un = dMA_span.array()/sn.array();
+
+ /* Calculate velocities on control points in the x,y,z directions
+ On each control point we project the velocities in the chordwise,
+ spanwise and normal directions to the x,y,z axes */
+ Eigen::MatrixXd uc(2*pan->nPanChord, pan->nPanSpan);
+ Eigen::MatrixXd vc(2*pan->nPanChord, pan->nPanSpan);
+ Eigen::MatrixXd wc(2*pan->nPanChord, pan->nPanSpan);
+
+ Eigen::Matrix3d Rmat_sys;
+ Eigen::Vector3d RHS_sys;
+ Eigen::VectorXd dummy_sys;
+
+ for (size_t i=0;i<2*pan->nPanChord; ++i)
+ {
+ for (size_t j=0; j<pan->nPanSpan; ++j)
+ {
+ Rmat_sys(0,0) = tmx(i,j);
+ Rmat_sys(0,1) = tmy(i,j);
+ Rmat_sys(0,2) = tmz(i,j);
+
+ Rmat_sys(1,0) = tnx(i,j);
+ Rmat_sys(1,1) = tny(i,j);
+ Rmat_sys(1,2) = tnz(i,j);
+
+ Rmat_sys(2,0) = pan->nx(i,j);
+ Rmat_sys(2,1) = pan->ny(i,j);
+ Rmat_sys(2,2) = pan->nz(i,j);
+
+ RHS_sys(0) = um(i,j);
+ RHS_sys(1) = un(i,j);
+ RHS_sys(2) = SA(i,j);
+
+ dummy_sys = Rmat_sys.partialPivLu().solve(RHS_sys);
+ uc(i,j)=dummy_sys(0);
+ vc(i,j)=dummy_sys(1);
+ wc(i,j)=dummy_sys(2);
+ }
+ }
+
+ uc = Eigen::MatrixXd::Constant(uc.rows(), uc.cols(), fs_velocity(0)) + uc;
+ vc = Eigen::MatrixXd::Constant(vc.rows(), vc.cols(), fs_velocity(1)) + vc;
+ wc = Eigen::MatrixXd::Constant(wc.rows(), wc.cols(), fs_velocity(2)) + wc;
+
+ /* Calculate presssure coefficient */
+ Eigen::MatrixXd dPhidt;
+ if (it==1){
+ Eigen::Map<Eigen::VectorXd> Phi1Flatten(Phi[1].data(),Phi[1].size());
+ Eigen::Map<Eigen::VectorXd> Phi0Flatten(Phi0.data(),Phi0.size());
+ dPhidt = (Phi1Flatten-Phi0Flatten)/dt;
+ }
+ else{
+ Eigen::Map<Eigen::VectorXd> Phi1Flatten(Phi[it].data(),Phi[it].size());
+ Eigen::Map<Eigen::VectorXd> Phi0Flatten(Phi[it-1].data(),Phi[it-1].size());
+ dPhidt = (Phi1Flatten - Phi0Flatten)/dt;
+ }
+
+ Eigen::MatrixXd cp_dum = (uc.array()*uc.array() + vc.array()*vc.array() + wc.array()*wc.array())*1/fs_airspeed/fs_airspeed - 2*dPhidt.array()*1/fs_airspeed/fs_airspeed;
+ cp = Eigen::MatrixXd::Constant(uc.rows(), uc.cols(), 1) - cp_dum;
+
+
+ /* Forces on panels */
+
+ std::cout << "cp" << std::endl;
+ std::cout << cp << std::endl;
+ std::cout << "pan->s" << std::endl;
+ std::cout << pan->s << std::endl;
+ std::cout << "pan->nx" << std::endl;
+ std::cout << pan->nx << std::endl;
+ std::cout << "pan->Surface" << std::endl;
+ std::cout << pan->Surface << std::endl;
+
+ Eigen::MatrixXd fx = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
+ Eigen::MatrixXd fy = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
+ Eigen::MatrixXd fz = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
+
+ std::cout << fx << std::endl;
+ std::cout << fy << std::endl;
+ std::cout << fz << std::endl;
+
+ double Fx = fx.sum();
+ //double Fy = fy.sum();
+ double Fz = fz.sum();
+
+ if (R>0){
+ CL(it) = Fz;
+ CD(it) = Fx;
+ }
+ else{
+ CL(it) = Fz-Fx*fs_alpha;
+ CD(it) = Fz*fs_alpha+Fx;
+ }
+ std::cout << " CL = " << CL(it) << std::endl;
+ it++;
+
+ }
+}
+
+
+ // int nChords = 100;
+ // double dt = nChords * pan->Getc0() / fs_velocity[0];
+
+ // /* Source strength */
+
+ // Eigen::MatrixXd SA = -(fs_velocity[0]*pan->nx/fs_beta + fs_velocity[1]*pan->ny + fs_velocity[2]*pan->nz);
+ // Eigen::MatrixXd dum = SA.transpose();
+ // Eigen::MatrixXd sigma = dum.reshaped(2*pan->nPanSpan*pan->nPanChord,1);
+
+ // pan->CreateWake(fs_velocity[0], fs_velocity[1], fs_velocity[2], fs_beta, dt);
+ // std::cout << "Wake created." << std::endl;
+
+ // pan->ComputeInfluenceCoeffs_wing();
+ // pan->ComputeInfluenceCoeffs_wake();
+ // std::cout << "Computed influence coefficients." << std::endl;
+
+ // pan->SetKutta();
+ // std::cout << "Imposed Kutta condition." << std::endl;
+
+ // /* Set up Dirichlet boundary condition */
+
+ // /* Horizontal concatenation */
+ // Eigen::MatrixXd LHS_dum1(pan->BN.rows(), pan->BN.cols()+pan->CN.cols());
+ // LHS_dum1 << pan->BN, pan->CN;
+
+ // /* Vertical concatenation */
+ // Eigen::MatrixXd LHS(LHS_dum1.rows()+pan->Kutta.rows(), LHS_dum1.cols());
+ // LHS << LHS_dum1,
+ // pan->Kutta;
+
+ // Eigen::MatrixXd RHS_dum1;
+ // RHS_dum1 = -pan->AN*sigma;
+
+ // Eigen::MatrixXd RHS(RHS_dum1.rows()+Eigen::MatrixXd::Zero(pan->nPanSpan, 1).rows(), RHS_dum1.cols());
+ // RHS << RHS_dum1,
+ // Eigen::MatrixXd::Zero(pan->nPanSpan, 1);
+
+ // /* Solve for strength */
+ // Eigen::VectorXd Mu = LHS.partialPivLu().solve(RHS);
+
+ // Eigen::MatrixXd MA_dum1;
+ // MA_dum1 = Mu.reshaped(pan->nPanSpan, 2*pan->nPanChord+1).transpose();
+
+ // /* Wake doublet strengths */
+
+ // Eigen::RowVectorXd Mw = MA_dum1.row(MA_dum1.rows()-1); // Last row of MA_dum1
+
+ // Eigen::VectorXd Mu2(2*pan->nPanChord*pan->nPanSpan);
+
+ // for (size_t i=0; i<2*pan->nPanChord*pan->nPanSpan; ++i)
+ // {
+ // Mu2.row(i) = Mu.row(i);
+ // }
+
+ // Eigen::MatrixXd MA(2*pan->nPanChord, MA_dum1.cols());
+
+ // for (size_t i=0; i<2*pan->nPanChord; ++i)
+ // {
+ // MA.row(i) = MA_dum1.row(i);
+ // }
+
+ // /* Calculate derivative of doublet strength in the chordwise direction */
+ // /* dX */
+ // Eigen::MatrixXd dXc1_chord = (pan->Xc.row(1)-pan->Xc.row(0));
+ // Eigen::MatrixXd dXc2_chord = (pan->Xc.block(2,0,2*pan->nPanChord,pan->Xc.cols()) - pan->Xc.block(0,0,2*pan->nPanChord-2,pan->Xc.cols()));
+ // Eigen::MatrixXd dXc3_chord = (pan->Xc.row(2*pan->nPanChord-1)-pan->Xc.row(2*pan->nPanChord-2));
+
+ // /* Vertical concatenation */
+ // Eigen::MatrixXd dXc_chord(dXc1_chord.rows()+dXc2_chord.rows()+dXc3_chord.rows(), dXc1_chord.cols());
+ // dXc_chord << dXc1_chord,
+ // dXc2_chord,
+ // dXc3_chord;
+
+ // /* dY */
+ // Eigen::MatrixXd dYc1_chord= (pan->Yc.row(1)-pan->Yc.row(0));
+ // Eigen::MatrixXd dYc2_chord= (pan->Yc.block(2,0,2*pan->nPanChord,pan->Yc.cols()) - pan->Yc.block(0,0,2*pan->nPanChord-2,pan->Yc.cols()));
+ // Eigen::MatrixXd dYc3_chord = (pan->Yc.row(2*pan->nPanChord-1)-pan->Yc.row(2*pan->nPanChord-2));
+
+ // /* Vertical concatenation */
+ // Eigen::MatrixXd dYc_chord(dYc1_chord.rows()+dYc2_chord.rows()+dYc3_chord.rows(), dYc1_chord.cols());
+ // dYc_chord << dYc1_chord,
+ // dYc2_chord,
+ // dYc3_chord;
+
+ // /* dZ */
+ // Eigen::MatrixXd dZc1_chord = (pan->Zc.row(1)-pan->Zc.row(0));
+ // Eigen::MatrixXd dZc2_chord = (pan->Zc.block(2,0,2*pan->nPanChord,pan->Zc.cols()) - pan->Zc.block(0,0,2*pan->nPanChord-2,pan->Zc.cols()));
+ // Eigen::MatrixXd dZc3_chord = (pan->Zc.row(2*pan->nPanChord-1)-pan->Zc.row(2*pan->nPanChord-2));
+
+ // /* Vertical concatenation */
+ // Eigen::MatrixXd dZc_chord(dZc1_chord.rows()+dZc2_chord.rows()+dZc3_chord.rows(), dZc1_chord.cols());
+ // dZc_chord << dZc1_chord,
+ // dZc2_chord,
+ // dZc3_chord;
+
+ // /* dMA_chord */
+ // Eigen::MatrixXd dMA1_chord = (MA.row(1)-MA.row(0));
+ // Eigen::MatrixXd dMA2_chord = (MA.block(2,0,2*pan->nPanChord,MA.cols()) - MA.block(0,0,2*pan->nPanChord-2,MA.cols()));
+ // Eigen::MatrixXd dMA3_chord = (MA.row(2*pan->nPanChord-1)-MA.row(2*pan->nPanChord-2));
+
+ // /* Vertical concatenation */
+ // Eigen::MatrixXd dMA_chord(dMA1_chord.rows()+dMA2_chord.rows()+dMA3_chord.rows(), dMA1_chord.cols());
+ // dMA_chord << dMA1_chord,
+ // dMA2_chord,
+ // dMA3_chord;
+
+ // /* Calculate chordwise tangent vectors */
+ // /* .array() member of MatrixXd allows the structure to be treated as ArrayXd for the operation.
+ // This allows product element by element instead of the classic matrix product of MatrixXd (aliasing is present for ArrayXd*ArrayXd)*/
+ // Eigen::MatrixXd sc_chord = (dXc_chord.array()*dXc_chord.array()+dYc_chord.array()*dYc_chord.array()+dZc_chord.array()*dZc_chord.array()).sqrt();
+
+ // Eigen::MatrixXd tmx_chord = dXc_chord.array()/sc_chord.array();
+ // Eigen::MatrixXd tmy_chord = dYc_chord.array()/sc_chord.array();
+ // Eigen::MatrixXd tmz_chord = dZc_chord.array()/sc_chord.array();
+ // /* Calculate the velocity in the chordwise direction */
+ // Eigen::MatrixXd um = dMA_chord.array()/sc_chord.array();
+
+
+ // /* Calculate derivative of doublet strength in the spanwise direction */
+ // /* dX */
+ // Eigen::MatrixXd dXc1_span = (pan->Xc.col(1)-pan->Xc.col(0));
+ // Eigen::MatrixXd dXc2_span = (pan->Xc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Xc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ // Eigen::MatrixXd dXc3_span = (pan->Xc.col(pan->nPanSpan-1)-pan->Xc.col(pan->nPanSpan-2));
+
+ // /* Horizontal concatenation */
+ // Eigen::MatrixXd dXc_span(dXc1_span.rows(), dXc1_span.cols()+dXc2_span.cols()+dXc3_span.cols());
+ // dXc_span << dXc1_span, dXc2_span, dXc3_span;
+
+ // /* dY */
+ // Eigen::MatrixXd dYc1_span = (pan->Yc.col(1)-pan->Yc.col(0));
+ // Eigen::MatrixXd dYc2_span = (pan->Yc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Yc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ // Eigen::MatrixXd dYc3_span = (pan->Yc.col(pan->nPanSpan-1)-pan->Yc.col(pan->nPanSpan-2));
+
+ // /* Horizontal concatenation */
+ // Eigen::MatrixXd dYc_span(dYc1_span.rows(), dYc1_span.cols()+dYc2_span.cols()+dYc3_span.cols());
+ // dYc_span << dYc1_span, dYc2_span, dYc3_span;
+
+ // /* dZ */
+ // Eigen::MatrixXd dZc1_span = (pan->Zc.col(1)-pan->Zc.col(0));
+ // Eigen::MatrixXd dZc2_span = (pan->Zc.block(0,2,2*pan->nPanChord,pan->nPanSpan) - pan->Zc.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ // Eigen::MatrixXd dZc3_span = (pan->Zc.col(pan->nPanSpan-1)-pan->Zc.col(pan->nPanSpan-2));
+
+ // /* Horizontal concatenation */
+ // Eigen::MatrixXd dZc_span(dZc1_span.rows(), dZc1_span.cols()+dZc2_span.cols()+dZc3_span.cols());
+ // dZc_span << dZc1_span, dZc2_span, dZc3_span;
+
+ // /* dMA_span */
+ // Eigen::MatrixXd dMA1_span = (MA.col(1)-MA.col(0));
+ // Eigen::MatrixXd dMA2_span = (MA.block(0,2,2*pan->nPanChord,pan->nPanSpan) - MA.block(0,0,2*pan->nPanChord,pan->nPanSpan-2));
+ // Eigen::MatrixXd dMA3_span = (MA.col(pan->nPanSpan-1)-MA.col(pan->nPanSpan-2));
+
+ // /* Horizontal concatenation */
+ // Eigen::MatrixXd dMA_span(dMA1_span.rows(), dMA1_span.cols()+dMA2_span.cols()+dMA3_span.cols());
+ // dMA_span << dMA1_span, dMA2_span, dMA3_span;
+
+ // /* Calculate spanwise tangent vectors */
+ // /* .array() member of MatrixXd allows the structure to be treated as ArrayXd for the operation.
+ // This allows product element by element instead of the classic matrix product of MatrixXd (aliasing is present for ArrayXd*ArrayXd)*/
+ // Eigen::MatrixXd sc_span = (dXc_span.array()*dXc_span.array()+dYc_span.array()*dYc_span.array()+dZc_span.array()*dZc_span.array()).sqrt();
+
+ // Eigen::MatrixXd tnx_span = dXc_span.array()/sc_span.array();
+ // Eigen::MatrixXd tny_span = dYc_span.array()/sc_span.array();
+ // Eigen::MatrixXd tnz_span = dZc_span.array()/sc_span.array();
+ // /* Calculate the velocity in the chordwise direction */
+ // Eigen::MatrixXd un = dMA_span.array()/sc_span.array();
+
+ // /* Calculate velocities on control points in the x,y,z directions
+ // On each control point we project the velocities in the chordwise,
+ // spanwise and normal directions to the x,y,z axes */
+ // Eigen::MatrixXd uc(2*pan->nPanChord, pan->nPanSpan);
+ // Eigen::MatrixXd vc(2*pan->nPanChord, pan->nPanSpan);
+ // Eigen::MatrixXd wc(2*pan->nPanChord, pan->nPanSpan);
+
+ // Eigen::Matrix3d Rmat_sys;
+ // Eigen::Vector3d RHS_sys;
+ // Eigen::VectorXd dummy_sys;
+
+ // for (size_t i=0;i<2*pan->nPanChord; ++i)
+ // {
+ // for (size_t j=0; j<pan->nPanSpan; ++j)
+ // {
+ // Rmat_sys(0,0) = tmx_chord(i,j);
+ // Rmat_sys(0,1) = tmy_chord(i,j);
+ // Rmat_sys(0,2) = tmz_chord(i,j);
+
+ // Rmat_sys(1,0) = tnx_span(i,j);
+ // Rmat_sys(1,1) = tny_span(i,j);
+ // Rmat_sys(1,2) = tnz_span(i,j);
+
+ // Rmat_sys(2,0) = pan->nx(i,j);
+ // Rmat_sys(2,1) = pan->ny(i,j);
+ // Rmat_sys(2,2) = pan->nz(i,j);
+
+ // RHS_sys(0) = um(i,j);
+ // RHS_sys(1) = un(i,j);
+ // RHS_sys(2) = SA(i,j); // or SA0?
+
+ // dummy_sys = Rmat_sys.partialPivLu().solve(RHS_sys);
+ // uc(i,j)=dummy_sys(0);
+ // vc(i,j)=dummy_sys(1);
+ // wc(i,j)=dummy_sys(2);
+ // }
+ // }
+
+ // /* Scale velocity to compressible space */
+ // uc/=fs_beta;
+ // uc+=Eigen::MatrixXd::Constant(uc.rows(),uc.cols(),fs_velocity(0));
+ // vc+=Eigen::MatrixXd::Constant(vc.rows(),vc.cols(),fs_velocity(1));
+ // wc+=Eigen::MatrixXd::Constant(wc.rows(),wc.cols(),fs_velocity(2));
+
+ // std::cout << " uc " << uc << std::endl;
+ // std::cout << " vc " << vc << std::endl;
+ // std::cout << " wc " << wc << std::endl;
+
+
+ // /* Calculate pressure coefficient using second order relationship */
+ // Eigen::MatrixXd bracket = Eigen::MatrixXd::Constant(uc.rows(),uc.cols(),1/2).array()*(uc.array()*uc.array()+vc.array()*vc.array()+wc.array()*wc.array()) - Eigen::MatrixXd::Constant(uc.rows(),uc.cols(),1/2*fs_airspeed*fs_airspeed).array();
+
+ // std::cout << "brack " << bracket << std::endl;
+
+ // Eigen::MatrixXd cp=(-bracket.array() + Eigen::MatrixXd::Constant(bracket.rows(), bracket.cols(),1/2/(fs_asound*fs_asound)).array() * bracket.array() * bracket.array()) * Eigen::MatrixXd::Constant(bracket.rows(), bracket.cols(),2/(fs_airspeed*fs_airspeed)).array();
+ // std::cout << " cp " << cp << std::endl;
+
+ // % Calculate total velocity
+ // Vtot2=uc.^2+vc.^2+wc.^2;
+
+ // % Calculate maximum velocity (giving Mach=1);
+ // Vmax=Q*sqrt(1+2/((gama-1)*Mach^2));
+ // % Calculate local speed of sound
+ // asound_local=asound*sqrt(1-(gama-1)/2*Mach^2*(Vtot2/Q^2-1));
+ // % Calculate local Mach number
+ // Mach_local=sqrt(Vtot2)./asound_local;
+ // % Calculate critical velocity (where local Mach number = 1)
+ // Vcrit=sqrt((gama-1)/(gama+1))*Vmax;
+ // iko=find(sqrt(Vtot2) >= Vcrit);
+ // if ~isempty(iko)
+ // disp('Warning! Sonic Mach number exceeded')
+ // super=1;
+ // end
+ // supervec(ialpha)=super;
+
+ // % Calculate forces on panels
+ // fx=-cp.*s0.*nx0/S;
+ // fy=-cp.*s0.*ny0/S;
+ // fz=-cp.*s0.*nz0/S;
+
+ // % Calculate total forces on wing in the x,y,z directions
+ // Fx=sum(sum(fx));
+ // Fy=sum(sum(fy));
+ // Fy_half=sum(sum(fy(:,n/2+1:n)));
+ // Fz=sum(sum(fz));
+
+ // % Calculate lift and drag coefficients
+ // CL(ialpha)=Fz*cos(alpha)-Fx*sin(alpha);
+ // CD(ialpha)=Fz*sin(alpha)+Fx*cos(alpha);
+ // CY(ialpha)=Fy;
+ // CY_half(ialpha)=Fy_half;
+
+ // % Calculate lift and drag distributions
+ // dy=Yp(1,2:n+1)-Yp(1,1:n);
+ // cc=(cvar(2:end)+cvar(1:end-1))/2;
+ // cfx=S*sum(fx,1)./dy./cc;
+ // cfy=S*sum(fy,1)./dy./cc;
+ // cfz=S*sum(fz,1)./dy./cc;
+ // Cl=cfz*cos(alpha)-cfx*sin(alpha);
+ // Cd=cfx*cos(alpha)+cfz*sin(alpha);
\ No newline at end of file
diff --git a/sdpm/src/pSDSolver.h b/sdpm/src/pSDSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..9180461797d056018252e3c4ca140d49a315f457
--- /dev/null
+++ b/sdpm/src/pSDSolver.h
@@ -0,0 +1,35 @@
+#ifndef PSDSOLVER_H
+#define PSDSOLVER_H
+
+#include "sdpm.h"
+#include "pPanels.h"
+#include <vector>
+#include <string>
+
+namespace sdpm
+{
+
+class DART_API SDSolver
+{
+
+private:
+
+
+
+public:
+
+Panels *pan;
+double fs_alpha, fs_Mach, fs_airspeed;
+unsigned long maxIt;
+Eigen::VectorXd CL;
+Eigen::VectorXd CD;
+Eigen::VectorXd CY;
+Eigen::MatrixXd cp;
+
+SDSolver(Panels *_pan, double _alpha, double _airspeed, unsigned long _maxIt);
+void Run();
+
+
+};
+} // namespace sdpm
+#endif //WSDSolver_H
diff --git a/sdpm/src/sdpm.h b/sdpm/src/sdpm.h
new file mode 100644
index 0000000000000000000000000000000000000000..c5dc54e405f2404dfcccd3a3b39822fc793d8822
--- /dev/null
+++ b/sdpm/src/sdpm.h
@@ -0,0 +1,44 @@
+/*
+ * 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.
+ */
+
+// global header of the "sdpm" module
+
+#ifndef sdpm_H
+#define sdpm_H
+
+#if defined(WIN32)
+#ifdef sdpm_EXPORTS
+#define DART_API __declspec(dllexport)
+#else
+#define DART_API __declspec(dllimport)
+#endif
+#else
+#define DART_API
+#endif
+
+#include "tbox.h"
+
+/**
+ * @brief Namespace for dart module
+ */
+namespace sdpm
+{
+// Viscous
+class Panels;
+class SDSolver;
+} // namespace sdpm
+
+#endif //SDPM_H
diff --git a/sdpm/tests/panels3d.py b/sdpm/tests/panels3d.py
new file mode 100644
index 0000000000000000000000000000000000000000..04bb3b13f3538e979c1c7987b1d0819fe1ab58a3
--- /dev/null
+++ b/sdpm/tests/panels3d.py
@@ -0,0 +1,27 @@
+import sdpm
+import math
+from matplotlib import pyplot as plt
+
+def main():
+
+ alpha = 2.*math.pi/180
+
+ nPanels = 2
+ airspeed = 50
+ maxIt = 10
+
+ wing = sdpm.Panels(nPanels)
+ wing.CreateWing()
+
+ solver = sdpm.SDSolver(wing, alpha, airspeed, maxIt)
+ print('Computation started (pyCheck)')
+
+ solver.Run()
+
+ plt.plot(solver.CL)
+ plt.show()
+
+ print('Computation ended (pyCheck)')
+
+if __name__ == "__main__":
+ main()
\ No newline at end of file
diff --git a/vii/_src/viiw.h b/vii/_src/viiw.h
index 806068743093274c2a5ba04e883db833fca3158f..0779519f043bd4dbbd87e725d764a15986524687 100644
--- a/vii/_src/viiw.h
+++ b/vii/_src/viiw.h
@@ -20,5 +20,4 @@
#include "wViscFluxes.h"
#include "wClosures.h"
#include "wViscMesh.h"
-#include "wInvBnd.h"
-#include "wCoupler.h"
\ No newline at end of file
+#include "wInvBnd.h"
\ No newline at end of file
diff --git a/vii/_src/viiw.i b/vii/_src/viiw.i
index 2d788b18118ef354402b03998f6ca3fc4181d2ca..4afee6e1640b86263f53cbc72fbd7a245edcbab1 100644
--- a/vii/_src/viiw.i
+++ b/vii/_src/viiw.i
@@ -49,5 +49,4 @@ threads="1"
%include "wViscFluxes.h"
%include "wClosures.h"
%include "wViscMesh.h"
-%include "wInvBnd.h"
-%include "wCoupler.h"
\ No newline at end of file
+%include "wInvBnd.h"
\ No newline at end of file
diff --git a/vii/pyVII/DAirfoil.py b/vii/pyVII/DAirfoil.py
index e8986c86036c040ebac6860fbd7b67afbeacd81c..70605b55cedbb79d6a893214bd451cb3341ff815 100755
--- a/vii/pyVII/DAirfoil.py
+++ b/vii/pyVII/DAirfoil.py
@@ -135,6 +135,4 @@ class Airfoil(Boundary):
uData = data[0:np.argmax(data[:,0])+1]
lData = data[np.argmax(data[:,0])+1:None]
lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
- print(np.shape(connectListNodes))
- quit()
return connectListNodes, connectListElems,[uData, lData]
diff --git a/vii/pyVII/vCoupler.py b/vii/pyVII/vCoupler.py
index 78d78c9237cb2fe0fd4530c789cf9d240b5a925f..5378573b2824206f1bf556e2d6aee5aa751e8493 100644
--- a/vii/pyVII/vCoupler.py
+++ b/vii/pyVII/vCoupler.py
@@ -193,8 +193,6 @@ class Coupler:
self.fMeshDict, _, _ = GroupConstructor().mergeVectors(dataIn, 0, 1)
fMeshDict = self.fMeshDict.copy()
if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
-
- print('STAG POINT MOVED')
# Initialize mesh
diff --git a/vii/src/wBLRegion.cpp b/vii/src/wBLRegion.cpp
index 0f83eccf4342444ca174461ca471b6967f776a3f..deb501f57579cd604c7e1b73e956b16beaee72d6 100644
--- a/vii/src/wBLRegion.cpp
+++ b/vii/src/wBLRegion.cpp
@@ -2,12 +2,8 @@
#include "wClosures.h"
#include "wViscMesh.h"
#include "wInvBnd.h"
-#include <iostream>
-#include <vector>
#include <cmath>
-
#include <iomanip>
-#include <deque>
using namespace vii;
@@ -35,6 +31,9 @@ BLRegion::~BLRegion()
// std::cout << "~BLRegion()\n";
} // end ~BLRegion
+/**
+ * @brief Set mesh in the global frame of reference
+ */
void BLRegion::SetMesh(std::vector<double> x,
std::vector<double> y,
std::vector<double> z)
@@ -63,6 +62,9 @@ void BLRegion::SetMesh(std::vector<double> x,
}
} // End SetMesh
+/**
+ * @brief Impose boundary condition at the stagnation point
+ */
void BLRegion::SetStagBC(double Re)
{
double dv0 = (invBnd->GetVt(1) - invBnd->GetVt(0)) / (mesh->GetLoc(1) - mesh->GetLoc(0));
@@ -105,6 +107,9 @@ void BLRegion::SetStagBC(double Re)
} // End SetStagBC
+/**
+ * @brief Impose boundary condition at the first wake point
+ */
void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
{
if (name!="wake")
@@ -145,6 +150,9 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
}
}
+/**
+ * @brief Compute blowing velocities at each station
+ */
void BLRegion::ComputeBlwVel()
{
DeltaStar[0] = U[0] * U[1];
@@ -157,7 +165,9 @@ void BLRegion::ComputeBlwVel()
}
}
-
+/**
+ * @brief Print solution at a given station
+ */
void BLRegion::PrintSolution(size_t iPoint) const
{
std::cout << "Pt " << iPoint << "Reg " << Regime[iPoint] << std::endl;
diff --git a/vii/src/wBLRegion.h b/vii/src/wBLRegion.h
index b109fdf7ce271c48f47458e8f1199b5f648b7a6d..e1b197966f08d53dca8b32500a931c746e8ca9cb 100644
--- a/vii/src/wBLRegion.h
+++ b/vii/src/wBLRegion.h
@@ -5,8 +5,6 @@
#include "wClosures.h"
#include "wViscMesh.h"
#include "wInvBnd.h"
-#include <vector>
-#include <string>
namespace vii
{
@@ -15,44 +13,44 @@ class DART_API BLRegion
{
private:
- unsigned int const nVar = 5;
+ unsigned int const nVar = 5; /* Number of variables of the partial differential problem */
/* Transition related variables */
- double nCrit = 9.0;
+ double nCrit = 9.0; /* Critical amplification factor */
public:
- ViscMesh *mesh;
- InvBnd *invBnd;
- Closures *closSolver;
+ ViscMesh *mesh; /* Boundary layer mesh */
+ InvBnd *invBnd; /* Inviscid quantities at the edge of the boundary layer */
+ Closures *closSolver; /* Closure relatrions module */
- std::string name; /* Name of the region */
+ std::string name; /* Name of the region */
/* Boundary layer variables */
- std::vector<double> U; /* Solution vector (theta, H, N, ue, Ct) */
- std::vector<double> Ret; /* Reynolds number based on the momentum thickness (theta) */
- std::vector<double> Cd; /* Local dissipation coefficient */
- std::vector<double> Cf; /* Local friction coefficient */
- std::vector<double> Hstar;
- std::vector<double> Hstar2;
- std::vector<double> Hk;
- std::vector<double> Cteq;
- std::vector<double> Us;
- std::vector<double> Delta;
+ std::vector<double> U; /* Solution vector (theta, H, N, ue, Ct) */
+ std::vector<double> Ret; /* Reynolds number based on the momentum thickness (theta) */
+ std::vector<double> Cd; /* Local dissipation coefficient */
+ std::vector<double> Cf; /* Local friction coefficient */
+ std::vector<double> Hstar; /* Kinetic energy shape parameter (thetaStar/theta) */
+ std::vector<double> Hstar2; /* Density shape parameter (deltaStarStar/theta) */
+ std::vector<double> Hk; /* Kinematic shape parameter (int(1-u/ue deta)) */
+ std::vector<double> Cteq; /* Equilibrium shear stress coefficient (turbulent BL) */
+ std::vector<double> Us; /* Equivalent normalized wall slip velocity */
+ std::vector<double> Delta; /* Boundary layer thickness */
- std::vector<double> DeltaStar;
+ std::vector<double> DeltaStar; /* Dispacement thickness (int(1-rho*u/rhoe*ue deta)) */
/* Transition related variables */
- std::vector<int> Regime;
- std::vector<double> Blowing;
+ std::vector<int> Regime; /* Laminar (0) or turbulent (1) regime */
+ std::vector<double> Blowing; /* Blowing velocity */
- unsigned int transMarker;
- double xtrF;
- double xtr;
+ unsigned int transMarker; /* Marker of the transition location */
+ double xtrF; /* Forced transition location */
+ double xtr; /* Transition location (=xtrF if forced) */
BLRegion(double _xtrF = -1.0);
~BLRegion();
diff --git a/vii/src/wClosures.cpp b/vii/src/wClosures.cpp
index 10b2349952d107b8005be1a4273453f8f58a9a70..49c1aa85bc2c1afb541e995f94a2c26fb4adf61d 100644
--- a/vii/src/wClosures.cpp
+++ b/vii/src/wClosures.cpp
@@ -1,19 +1,17 @@
#include "wClosures.h"
+#include <cmath>
#include <iostream>
-#include <math.h>
-#include <string>
-using namespace tbox;
using namespace vii;
// This is an empty constructor
Closures::Closures() {}
-Closures::~Closures()
-{
- // std::cout << "~Closures()\n";
-} // end ~Closures
+Closures::~Closures(){}
+/**
+ * @brief Laminar closure relations
+ */
std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const
{
@@ -49,7 +47,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
{
Hstar = 1.528 + 0.015*Hks*Hks/Hk;
}
- // Whitfield's minor additional compressibility correction
+ /* Whitfield's minor additional compressibility correction */
Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me);
/* Friction coefficient */
@@ -106,6 +104,9 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
return ClosureVars;
}
+/**
+ * @brief Turbulent closure relations
+ */
std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const
{
/*if (std::isnan(Ret))
diff --git a/vii/src/wDartCoupler.cpp b/vii/src/wDartCoupler.cpp
deleted file mode 100644
index cd40af47155f42ede0f8215ba6116422f3968025..0000000000000000000000000000000000000000
--- a/vii/src/wDartCoupler.cpp
+++ /dev/null
@@ -1,172 +0,0 @@
-#include "wDartCoupler.h"
-#include "../../dart/src/wNewton.h"
-#include "../../dart/src/wBlowing.h"
-#include "wViscSolver.h"
-
-using namespace vii;
-
-DartCoupler::DartCoupler(const Newton *_invSolver, ViscSolver *_viscSolver, const std::vector<*Blowing> _blw, unsigned long _nCouplIter = 150, double _tol = 1e-4)
-{
- invSolver = _invSolver;
- viscSolver = _viscSolver;
- blw = _blw;
-
- maxCouplIter = _nCouplIter;
- tol = _tol;
-
- Initialize();
-}
-DartCoupler::~DartCoupler()
-
-void DartCoupler::Run(){
- double cdPrev = 0.;
- double error = 1.;
-
- unsigned long nIters;
-
- while (nIters < maxCouplIter){
-
- invSolver.run();
-
- ImposeInviscidBC();
-
- viscSolver.Run();
-
- error = (viscSolver.Cdt - cdPrev) / viscSolver.Cdt;
-
- if (error < tol){
- break;
- }
- else{
- continue;
- }
-
- cdPrev = viscSolver.Cdt;
-
- ImposeBlowing();
-
- nIters++;
-
- }
-
- else{
- std::cout << "Maximum number of " << maxCouplIter << " reached." << std::endl;
- }
-}
-void DartCoupler::ImposeInviscidBC(){
-
- for (int iReg=0; iReg < blw.size(); ++iReg){
-
- for (auto iNode : blw[iReg]->nodes){
- invSolver->U(iNode->row)
- }
- }
-}
-
-void DartCoupler::ImposeBlowing()
-
-void DartCoupler::Initialize(){
-
- /* Create connectivity list */
-
- std::vector<std::vector<unsigned int>,2> connectListNodes;
- std::vector<std::vector<unsigned int>,2> connectListElems;
-
- for (size_t iReg=0; iReg < blw.size(); ++iReg)
- {
- connectListNodes(iReg).resize(blw(iReg)->boundary->nodes.size(), 0);
- connectListElems(iReg).resize(blw(iReg)->boudnary->tag->elems.size(), 0);
- }
-
-
- /* Find stagnation point */
- int idxStag =
-
-
- N1 = np.zeros(self.nN, dtype=int) # Node number
- connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
- connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
- data = np.zeros((self.boundary.nodes.size(), 10))
- i = 0
- for n in self.boundary.nodes:
- data[i,0] = n.no
- data[i,1] = n.pos[0]
- data[i,2] = n.pos[1]
- data[i,3] = n.pos[2]
- data[i,4] = self.v[i,0]
- data[i,5] = self.v[i,1]
- data[i,6] = self.Me[i]
- data[i,7] = self.rhoe[i]
- data[i,8] = self.deltaStar[i]
- data[i,9] = self.xx[i]
- i += 1
- # Table containing the element and its nodes
- eData = np.zeros((self.nE,3), dtype=int)
- for i in range(0, self.nE):
- eData[i,0] = self.boundary.tag.elems[i].no
- eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
- eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
- # Find the stagnation point
- if couplIter == 0:
- self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
- idxStag = self.idxStag
- else:
- idxStag = self.idxStag
- #idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
- globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
- # Find TE nodes (assuming suction side element is above pressure side element in the y direction)
- idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
- idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
- idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
- ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
- ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
- if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
- upperGlobTE = int(data[idxTE[0],0])
- lowerGlobTE = int(data[idxTE[1],0])
- else:
- upperGlobTE = int(data[idxTE[1],0])
- lowerGlobTE = int(data[idxTE[0],0])
- # Connectivity
- connectListElems[0] = np.where(eData[:,1] == globStag)[0]
- N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
- connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
- i = 1
- upperTE = 0
- lowerTE = 0
- # Sort the suction part
- while upperTE == 0:
- N1[i] = eData[connectListElems[i-1],2] # Second node of the element
- connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
- connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
- if eData[connectListElems[i],2] == upperGlobTE:
- upperTE = 1
- i += 1
- # Sort the pressure side
- connectListElems[i] = np.where(eData[:,2] == globStag)[0]
- connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
- N1[i] = eData[connectListElems[i],2]
- while lowerTE == 0:
- N1[i+1] = eData[connectListElems[i],1] # First node of the element
- connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
- connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
- if eData[connectListElems[i+1],1] == lowerGlobTE:
- lowerTE = 1
- i += 1
- connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
- data[:,0] = data[connectListNodes,0]
- data[:,1] = data[connectListNodes,1]
- data[:,2] = data[connectListNodes,2]
- data[:,3] = data[connectListNodes,3]
- data[:,4] = data[connectListNodes,4]
- data[:,5] = data[connectListNodes,5]
- data[:,6] = data[connectListNodes,6]
- data[:,7] = data[connectListNodes,7]
- data[:,8] = data[connectListNodes,8]
- data[:,9] = data[connectListNodes,9]
- # Separated upper and lower part
- data = np.delete(data,0,1)
- uData = data[0:np.argmax(data[:,0])+1]
- lData = data[np.argmax(data[:,0])+1:None]
- lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
- return connectListNodes, connectListElems,[uData, lData]
-}
diff --git a/vii/src/wDartCoupler.h b/vii/src/wDartCoupler.h
deleted file mode 100644
index 8910ecda39b957620d86509fe04c14d254ea6797..0000000000000000000000000000000000000000
--- a/vii/src/wDartCoupler.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef WDARTCOUPLER_H
-#define WDARTCOUPLER_H
-
-#include "vii.h"
-#include "../../dart/src/wNewton.h"
-#include "../../dart/src/wBlowing.h"
-#include "wViscSolver.h"
-
-namespace vii
-{
-
-class DART_API DartCoupler
-{
-
-private:
- unsigned long maxCouplIter;
- double tol;
-
-
-public:
-
- Newton *invSolver;
- ViscSolver *viscSolver;
- std::vector<Blowing> *blw;
-
-
- DartCoupler();
- ~DartCoupler();
- Run();
- ImposeBlowing();
- ImposeInviscidBC();
-
-};
-} // namespace dart
-#endif //WDARTCOUPLER_H
diff --git a/vii/src/wInvBnd.cpp b/vii/src/wInvBnd.cpp
index 8d8722792b13f37e4afb4a033f2e8da8d05a60af..7c28ab295d51809c854b98880090b35a98b14b38 100644
--- a/vii/src/wInvBnd.cpp
+++ b/vii/src/wInvBnd.cpp
@@ -11,13 +11,6 @@
/* --- Project includes -------------------------------------------------------------- */
#include "wInvBnd.h"
-#include <iostream>
-#include <vector>
-#include <unordered_set>
-#include <iomanip>
-#include <fstream>
-#include <iostream>
-#include <iomanip>
/* --- Namespaces -------------------------------------------------------------------- */
diff --git a/vii/src/wInvBnd.h b/vii/src/wInvBnd.h
index b3f714cad2bc913f878b01546a2a04c430e26817..6cfb14fbd8601221703b688b20031ad7c254feee 100644
--- a/vii/src/wInvBnd.h
+++ b/vii/src/wInvBnd.h
@@ -2,7 +2,6 @@
#define WINVBND_H
#include "vii.h"
-#include "wBLRegion.h"
#include <vector>
#include <iostream>
diff --git a/vii/src/wTimeSolver.cpp b/vii/src/wTimeSolver.cpp
index dd420136a76b9353f2b23017068625233f1c60e5..be6e0e7570c1ace5c50c68e9364ae9cd5bf05296 100644
--- a/vii/src/wTimeSolver.cpp
+++ b/vii/src/wTimeSolver.cpp
@@ -1,11 +1,7 @@
#include "wTimeSolver.h"
#include "wBLRegion.h"
#include "wViscFluxes.h"
-#include "wInvBnd.h"
#include <Eigen/Dense>
-#include <iostream>
-#include <vector>
-#include <cmath>
using namespace Eigen;
using namespace tbox;
diff --git a/vii/src/wTimeSolver.h b/vii/src/wTimeSolver.h
index af485d4cdb6686206da9a04252711821ba256569..e80e714dfcd3853aac9b0dee52db53448c2eb0c4 100644
--- a/vii/src/wTimeSolver.h
+++ b/vii/src/wTimeSolver.h
@@ -3,8 +3,6 @@
#include "vii.h"
#include "wViscFluxes.h"
-#include <vector>
-#include <string>
namespace vii
{
diff --git a/vii/src/wViscFluxes.cpp b/vii/src/wViscFluxes.cpp
index 7fe7947fbfab3dd2f0b0ff734502c6c0b927dce6..1b667ba5dfa3a88f3c768a753c0001c078623bcb 100644
--- a/vii/src/wViscFluxes.cpp
+++ b/vii/src/wViscFluxes.cpp
@@ -2,12 +2,7 @@
#include "wViscFluxes.h"
#include "wBLRegion.h"
-#include "wInvBnd.h"
-#include <iostream>
#include <Eigen/Dense>
-#include <cmath>
-#include <algorithm>
-#include <string>
using namespace vii;
using namespace Eigen;
diff --git a/vii/src/wViscFluxes.h b/vii/src/wViscFluxes.h
index a60bf35f34c0e2099b608e539694c2f4fd196a44..823a070e298f8c407d7d06b8af093a1424eb9988 100644
--- a/vii/src/wViscFluxes.h
+++ b/vii/src/wViscFluxes.h
@@ -4,8 +4,6 @@
#include "vii.h"
#include "wBLRegion.h"
#include <Eigen/Dense>
-#include <vector>
-#include <string>
namespace vii
{
diff --git a/vii/src/wViscMesh.cpp b/vii/src/wViscMesh.cpp
index 9cb054452bdf98b13c0f193d6493d7ab10d08d37..82f742274474f98069fbe55731186d0c5dcc8e64 100644
--- a/vii/src/wViscMesh.cpp
+++ b/vii/src/wViscMesh.cpp
@@ -11,8 +11,6 @@
/* --- Project includes -------------------------------------------------------------- */
#include "wViscMesh.h"
-#include <iostream>
-#include <vector>
#include <cmath>
/* --- Namespaces -------------------------------------------------------------------- */
diff --git a/vii/src/wViscMesh.h b/vii/src/wViscMesh.h
index 0531097bc4bcbec7c300c2278860b0f9109134b7..5dbec9743957e00eca017f9c7905769780a3ae19 100644
--- a/vii/src/wViscMesh.h
+++ b/vii/src/wViscMesh.h
@@ -3,8 +3,6 @@
#include "vii.h"
#include "wBLRegion.h"
-#include <vector>
-#include <iostream>
namespace vii
{
diff --git a/vii/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp
index e4270fed1af41494016b3686cacc5a1632cf4e54..c7dcb053c7c74a6f7f95b158e1e62178cc0080e5 100644
--- a/vii/src/wViscSolver.cpp
+++ b/vii/src/wViscSolver.cpp
@@ -13,10 +13,7 @@
#include "wViscSolver.h"
#include "wBLRegion.h"
#include "wTimeSolver.h"
-#include <iostream>
#include <iomanip>
-#include <vector>
-#include <cmath>
/* --- Namespaces -------------------------------------------------------------------- */
@@ -100,8 +97,8 @@ ViscSolver::~ViscSolver()
*/
int ViscSolver::Run(unsigned int const couplIter)
{
- bool lockTrans=false;
- int pointExitCode = 0;
+ bool lockTrans=false; /* Flagging activation of transition routines */
+ int pointExitCode = 0; /* Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter) */
solverExitCode = 0;
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
@@ -126,7 +123,7 @@ int ViscSolver::Run(unsigned int const couplIter)
Sections[iSec][0]->xtr = 1.; /* Upper side initial xtr */
Sections[iSec][1]->xtr = 1.; /* Lower side initial xtr */
- Sections[iSec][2]->xtr = 0.; /* Wake initial xtr (full turbulent) */
+ Sections[iSec][2]->xtr = 0.; /* Wake initial xtr (fully turbulent) */
Sections[iSec][0]->SetStagBC(Re);
Sections[iSec][1]->SetStagBC(Re);
@@ -137,7 +134,7 @@ int ViscSolver::Run(unsigned int const couplIter)
if (verbose>0){std::cout << Sections[iSec][iRegion]->name << " ";}
- /* Check if we need to enter safe mode */
+ /* Check if safe mode is required */
if (couplIter == 0 || Sections[iSec][iRegion]->mesh->GetDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
{
@@ -165,6 +162,7 @@ int ViscSolver::Run(unsigned int const couplIter)
Sections[iSec][iRegion]->invBnd->SetVtOld(ueOld);
}
+ /* Keyword annoncing iteration safety level */
if (verbose>1){std::cout << "restart. ";}
}
@@ -174,6 +172,7 @@ int ViscSolver::Run(unsigned int const couplIter)
tSolver->SetitFrozenJac(5);
tSolver->SetinitSln(false);
+ /* Keyword annoncing iteration safety level */
if (verbose>1){std::cout << "update. ";}
}
@@ -276,7 +275,7 @@ int ViscSolver::Run(unsigned int const couplIter)
}
/**
- * @brief Check whether instabilities are growing in the laminar boundary layer.
+ * @brief Check whether or not instabilities are growing in the laminar boundary layer.
*/
void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
{
@@ -302,14 +301,16 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
unsigned int nVar = bl->GetnVar();
- // Save laminar solution.
+ /* Save laminar solution. */
+
std::vector<double> lamSol(nVar, 0.0);
for (size_t k = 0; k < lamSol.size(); ++k)
{
lamSol[k] = bl->U[iPoint * nVar + k];
}
- // Set up turbulent portion boundary condition.
+ /* Set up turbulent portion boundary condition. */
+
bl->transMarker = iPoint; /* Save transition marker */
/* Loop over remaining points in the region */
@@ -319,7 +320,6 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
}
/* Compute transition location */
- //std::cout << "Trans " << bl->mesh->Getx(iPoint) << std::endl;
bl->xtr = (bl->GetnCrit() - bl->U[(iPoint - 1) * nVar + 2]) * ((bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1)) / (bl->U[iPoint * nVar + 2] - bl->U[(iPoint - 1) * nVar + 2])) + bl->mesh->Getx(iPoint - 1);
/* Percentage of the subsection corresponding to a turbulent flow */
diff --git a/vii/tests/bli2.py b/vii/tests/bli2.py
index c46a10a8e9a8fc058a3c827bd315b07aff9bf728..cf3b47694b2f44608f591eb48f8f30295baca732 100644
--- a/vii/tests/bli2.py
+++ b/vii/tests/bli2.py
@@ -57,7 +57,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5.*math.pi/180
+ alpha = 2.*math.pi/180
M_inf = 0.
CFL0 = 1