diff --git a/fpm/_src/fpmw.h b/fpm/_src/fpmw.h
index 59e1f6abbbd1fc34219e567d8808a5105443ffc5..0045215c39bf810698a597711cbbcee3968d4540 100644
--- a/fpm/_src/fpmw.h
+++ b/fpm/_src/fpmw.h
@@ -15,6 +15,8 @@
*/
#include "fBody.h"
+#include "fBuilder.h"
#include "fProblem.h"
+#include "fSolver.h"
#include "fTimers.h"
#include "fWake.h"
diff --git a/fpm/_src/fpmw.i b/fpm/_src/fpmw.i
index b68248d2b3c6045c12dbdb11b2cb1cf1a4b0b304..c09a59a21013e225161fc86d546121ac09cd20b2 100644
--- a/fpm/_src/fpmw.i
+++ b/fpm/_src/fpmw.i
@@ -74,5 +74,13 @@ def __getitem__(self, name):
%include "fWake.h"
%shared_ptr(fpm::Problem);
-%immutable fpm::Problem::msh; // avoids the creation of the setter method
+%immutable fpm::Problem::msh;
%include "fProblem.h"
+
+%shared_ptr(fpm::Builder);
+%immutable fpm::Builder::pbl;
+%include "fBuilder.h"
+
+%shared_ptr(fpm::Solver);
+%immutable fpm::Solver::aic;
+%include "fSolver.h"
diff --git a/fpm/src/CMakeLists.txt b/fpm/src/CMakeLists.txt
index 54ef33c32d192cb9fe7d103c895a8d3d76dec79d..c583cfd623cfc06babfe679eaa52d0d7fe66bcfe 100644
--- a/fpm/src/CMakeLists.txt
+++ b/fpm/src/CMakeLists.txt
@@ -21,12 +21,6 @@ MACRO_DebugPostfix(fpm)
TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${PROJECT_SOURCE_DIR}/fpm/src)
-# -- TBB
-FIND_PACKAGE(TBB REQUIRED)
-SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${TBB_CXX_FLAGS_DEBUG}")
-TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${TBB_INCLUDE_DIRS})
-TARGET_LINK_LIBRARIES(fpm ${TBB_LIBRARIES})
-
# -- Eigen
FIND_PACKAGE(EIGEN 3.3.4 REQUIRED)
TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${EIGEN_INCLUDE_DIRS})
diff --git a/fpm/src/fBody.cpp b/fpm/src/fBody.cpp
index 10e9200919399c52368a26c8e4866072772ca6c7..571fb5804c2c95b8c990a0611d1e87fcaf6791aa 100644
--- a/fpm/src/fBody.cpp
+++ b/fpm/src/fBody.cpp
@@ -30,6 +30,17 @@ using namespace fpm;
Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
std::vector<std::string> const &teIds, double xF) : Group(_msh, id), Cl(0), Cd(0), Cs(0), Cm(0)
{
+ // Get nodes
+ for (auto e : tag->elems)
+ for (auto n : e->nodes)
+ nodes.push_back(n);
+ std::sort(nodes.begin(), nodes.end());
+ nodes.erase(std::unique(nodes.begin(), nodes.end()), nodes.end());
+ // Associate each node to its adjacent elements
+ for (auto e : tag->elems)
+ for (auto n : e->nodes)
+ neMap[n].push_back(e);
+
// If wakes are requested, check if they are already available, otherwise create them
bool hasChanged = false;
size_t j = 0;
@@ -94,9 +105,9 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
}
// Size load vectors
- cLoadX.resize(tag->elems.size());
- cLoadY.resize(tag->elems.size());
- cLoadZ.resize(tag->elems.size());
+ cLoadX.resize(nodes.size());
+ cLoadY.resize(nodes.size());
+ cLoadZ.resize(nodes.size());
}
/**
diff --git a/fpm/src/fBody.h b/fpm/src/fBody.h
index 7cdea5ad3a2dd2b152f8ff2f60662d2cffeb44b0..181af38f670c9699a95aca280a52f37cec84b47c 100644
--- a/fpm/src/fBody.h
+++ b/fpm/src/fBody.h
@@ -20,25 +20,28 @@
#include "fpm.h"
#include "wGroup.h"
#include <vector>
+#include <map>
namespace fpm
{
/**
- * @brief Manage a body immersed in the fluid
+ * @brief Manage a body immersed in the fluid (lifting surface)
* @authors Adrien Crovato
*/
class FPM_API Body : public tbox::Group
{
public:
- double Cl; ///< lift coefficient
- double Cd; ///< drag coefficient
- double Cs; ///< sideforce coefficient
- double Cm; ///< pitch moment coefficient (positive nose-up)
- std::vector<double> cLoadX; ///< x-component of aerodynamic load coefficient on boundary
- std::vector<double> cLoadY; ///< y-component of aerodynamic load coefficient on boundary
- std::vector<double> cLoadZ; ///< z-component of aerodynamic load coefficient on boundary
- std::vector<Wake *> wakes; ///< wake(s) attached to the lifting surface
+ double Cl; ///< lift coefficient
+ double Cd; ///< drag coefficient
+ double Cs; ///< sideforce coefficient
+ double Cm; ///< pitch moment coefficient (positive nose-up)
+ std::vector<double> cLoadX; ///< x-component of aerodynamic load normalized by dynamic pressure
+ std::vector<double> cLoadY; ///< y-component of aerodynamic load normalized by dynamic pressure
+ std::vector<double> cLoadZ; ///< z-component of aerodynamic load normalized by dynamic pressure
+ std::vector<tbox::Node *> nodes; ///< nodes of the surface
+ std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< map between nodes and adjacent elements
+ std::vector<Wake *> wakes; ///< wake(s) attached to the lifting surface
Body(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::vector<std::string> const &teNames, double xF);
~Body() { std::cout << "~Body()\n"; }
diff --git a/fpm/src/fBuilder.cpp b/fpm/src/fBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc7bfddd435976f80a5a42f9d13986a7ea97eef5
--- /dev/null
+++ b/fpm/src/fBuilder.cpp
@@ -0,0 +1,162 @@
+/*
+ * 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 "fBuilder.h"
+#include "fProblem.h"
+//#include "fField.h"
+#include "fBody.h"
+#include "fWake.h"
+
+#include "wTag.h"
+#include "wElement.h"
+#include "wMem.h"
+
+#define PI 3.14159
+
+using namespace tbox;
+using namespace fpm;
+
+/**
+ * @brief Initialize the solver and perform sanity checks
+ * @authors Adrien Crovato
+ */
+Builder::Builder(std::shared_ptr<Problem> _pbl) : pbl(_pbl), valid(false)
+{
+ // Check the problem and update element memory
+ pbl->check();
+ pbl->updateMem();
+
+ // Setup AIC matrices
+ nF = 0; //pbl->field->tag->elems.size();
+ nP = 0;
+ for (auto body : pbl->bodies)
+ nP += body->tag->elems.size();
+ A.resize(nP, nP);
+ B.resize(nP, nP);
+ C.resize(nP, nF);
+ Cx.resize(nP, nF);
+ Cy.resize(nP, nF);
+ Cz.resize(nP, nF);
+}
+
+/**
+ * @brief Fill AIC matrices
+ */
+void Builder::run()
+{
+ // global->local id
+ std::map<Element *, int> rows;
+ int i = 0;
+ for (auto body : pbl->bodies)
+ for (auto e : body->tag->elems)
+ {
+ rows[e] = i;
+ i++;
+ }
+
+ // body to body
+ for (auto ibody : pbl->bodies)
+ for (auto ei : ibody->tag->elems)
+ {
+ // body panels
+ for (auto jbody : pbl->bodies)
+ for (auto ej : jbody->tag->elems)
+ {
+ A(rows.at(ei), rows.at(ej)) = mu(ei, ej);
+ B(rows.at(ei), rows.at(ej)) = tau(ei, ej);
+ }
+ // wake panels
+ for (auto jbody : pbl->bodies)
+ for (auto wake : jbody->wakes)
+ for (auto ew : wake->tag->elems)
+ {
+ double wmu = mu(ei, ew);
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).first)) += wmu;
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).second)) -= wmu;
+ }
+ }
+
+ // field to body
+ //for (auto body : pbl->bodies)
+ // for (auto e : body->tag->elems)
+ // for (size_t j = 0; j < field->tag->elems.size(); ++j)
+ // {
+ // C(rows.at(ei), j) = sigma(ei, field->tag->elems[j]);
+ // Cx(rows.at(ei), j) = sigmaX(ei, field->tag->elems[j]);
+ // Cy(rows.at(ei), j) = sigmaY(ei, field->tag->elems[j]);
+ // Cz(rows.at(ei), j) = sigmaY(ei, field->tag->elems[j]);
+ // }
+
+ valid = true;
+}
+
+/**
+ * @brief Compute doublet AIC from body panel ej to body panel ei
+ */
+double Builder::mu(Element *ei, Element *ej)
+{
+ if (ei == ej)
+ return 0.5;
+ else
+ return 0;
+}
+/**
+ * @brief Compute source AIC from body panel ej to body panel ei
+ */
+double Builder::tau(Element *ei, Element *ej)
+{
+ if (ei == ej)
+ return 0;
+ else
+ return 0;
+}
+/**
+ * @brief Compute source AIC from field panel ej to body panel ei
+ */
+double Builder::sigma(Element *ei, Element *ej)
+{
+ return 0;
+}
+/**
+ * @brief Compute source x-velocity AIC from field panel ej to body panel ei
+ */
+double Builder::sigmaX(Element *ei, Element *ej)
+{
+ return 0;
+}
+/**
+ * @brief Compute source y-velocity AIC from field panel ej to body panel ei
+ */
+double Builder::sigmaY(Element *ei, Element *ej)
+{
+ return 0;
+}
+/**
+ * @brief Compute source z-velocity AIC from field panel ej to body panel ei
+ */
+double Builder::sigmaZ(Element *ei, Element *ej)
+{
+ return 0;
+}
+
+void Builder::write(std::ostream &out) const
+{
+ out << "fpm::Builder with AIC matrices:"
+ << "\n\tA(" << A.rows() << "," << A.cols() << "), B(" << B.rows() << "," << B.cols() << ")"
+ << "\n\tC(" << C.rows() << "," << C.cols() << ")"
+ << "\n\tCx(" << Cx.rows() << "," << Cx.cols() << "), Cy(" << Cy.rows() << "," << Cy.cols() << "), Cz(" << Cz.rows() << "," << Cz.cols() << ")"
+ << "\nwith\t" << *pbl;
+}
diff --git a/fpm/src/fBuilder.h b/fpm/src/fBuilder.h
new file mode 100644
index 0000000000000000000000000000000000000000..bb290cc1eda4eb4ec983a34fea48bb60edd7b556
--- /dev/null
+++ b/fpm/src/fBuilder.h
@@ -0,0 +1,71 @@
+/*
+ * 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.
+ */
+
+#ifndef FBUILDER_H
+#define FBUILDER_H
+
+#include "fpm.h"
+#include "wObject.h"
+
+#include <iostream>
+#include <vector>
+#include <memory>
+#include <Eigen/Dense>
+
+namespace fpm
+{
+
+/**
+ * @brief Aerodynamic Influence Coefficients builder
+ * @authors Adrien Crovato
+ * @todo add symmetry support
+ * @todo implement AIC computation
+ */
+class FPM_API Builder : public fwk::wSharedObject
+{
+public:
+ std::shared_ptr<Problem> pbl; ///< problem definition
+
+ bool valid; ///< matrices are up-to-date
+ int nP; ///< number of body panels
+ int nF; ///< number of field panels
+ Eigen::MatrixXd A; ///< body-body doublet matrix
+ Eigen::MatrixXd B; ///< body-body source matrix
+ Eigen::MatrixXd C; ///< field-body source matrix
+ Eigen::MatrixXd Cx; ///< field-body x-velocity source matrix
+ Eigen::MatrixXd Cy; ///< field-body y-velocity source matrix
+ Eigen::MatrixXd Cz; ///< field-body z-velocity source matrix
+
+ Builder(std::shared_ptr<Problem> _pbl);
+ ~Builder() { std::cout << "~Builder()\n"; }
+
+ void run();
+
+#ifndef SWIG
+ virtual void write(std::ostream &out) const override;
+#endif
+
+private:
+ double mu(tbox::Element *ei, tbox::Element *ej);
+ double tau(tbox::Element *ei, tbox::Element *ej);
+ double sigma(tbox::Element *ei, tbox::Element *ej);
+ double sigmaX(tbox::Element *ei, tbox::Element *ej);
+ double sigmaY(tbox::Element *ei, tbox::Element *ej);
+ double sigmaZ(tbox::Element *ei, tbox::Element *ej);
+};
+
+} // namespace fpm
+#endif //FBUILDER_H
diff --git a/fpm/src/fProblem.cpp b/fpm/src/fProblem.cpp
index 5ac54291b181be9abdb8f9c0436e77302dfed82e..c2b666e9943bed48315d94c8f9a5fb54c3e46de2 100644
--- a/fpm/src/fProblem.cpp
+++ b/fpm/src/fProblem.cpp
@@ -19,12 +19,12 @@
#include "wTag.h"
#include "wElement.h"
#include "wMem.h"
+#include "wCache.h"
using namespace tbox;
using namespace fpm;
-Problem::Problem(std::shared_ptr<MshData> _smsh, double aoa, double aos, double sref, double cref,
- double xref, double yref, double zref) : smsh(_smsh), alpha(aoa), beta(aos),
- sR(sref), cR(cref)
+Problem::Problem(std::shared_ptr<MshData> _msh, double aoa, double aos, double mch,
+ double sref, double cref, double xref, double yref, double zref) : msh(_msh), alpha(aoa), beta(aos), mach(mch), sR(sref), cR(cref)
{
xR(0) = xref;
xR(1) = yref;
@@ -39,6 +39,62 @@ void Problem::add(std::shared_ptr<Body> b)
bodies.push_back(b);
}
+/**
+ * @brief Compute velocity at element
+ */
+Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &phi)
+{
+ Eigen::Vector3d v(0, 0, 0);
+ // @todo to be implemented if usefull
+ return v;
+ std::cout << std::endl;
+}
+
+/**
+ * @brief Compute density at element
+ */
+double Problem::Rho(tbox::Element &e, std::vector<double> const &phi)
+{
+ // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
+ Eigen::Vector3d u = U(e, phi);
+ double a = 1 + 0.2 * mach * mach * (1 - u.dot(u));
+ return sqrt(a * a * a * a * a);
+}
+
+/**
+ * @brief Compute Mach number at element
+ */
+double Problem::Mach(tbox::Element &e, std::vector<double> const &phi)
+{
+ if (mach == 0)
+ return 0;
+ else
+ {
+ // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
+ Eigen::Vector3d u = U(e, phi);
+ double a2 = 1 / (mach * mach) + 0.2 * (1 - u.dot(u)); // speed of sound squared
+ return u.norm() / sqrt(a2);
+ }
+}
+
+/**
+ * @brief Compute pressure coefficient at element
+ */
+double Problem::Cp(tbox::Element &e, std::vector<double> const &phi)
+{
+ Eigen::Vector3d u = U(e, phi);
+ if (mach == 0)
+ {
+ return 1 - u.dot(u);
+ }
+ else
+ {
+ //particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1);
+ double a = 1 + 0.2 * mach * mach * (1 - u.dot(u));
+ return 2 / (1.4 * mach * mach) * (sqrt(a * a * a * a * a * a * a) - 1);
+ }
+}
+
/**
* @brief Update the elements memory (Jacobian)
*/
diff --git a/fpm/src/fProblem.h b/fpm/src/fProblem.h
index 54cf0159924a8e100ad81a1ca01fc77b90669efb..7c8714537384f63befa5552b0dd923cc3e5270c9 100644
--- a/fpm/src/fProblem.h
+++ b/fpm/src/fProblem.h
@@ -28,26 +28,36 @@ namespace fpm
/**
* @brief Manage the problem
* @authors Adrien Crovato
+ * @todo check and implement velocity computation
*/
class FPM_API Problem : public fwk::wSharedObject
{
public:
- std::shared_ptr<tbox::MshData> smsh; ///< surface mesh data structure
+ std::shared_ptr<tbox::MshData> msh; ///< mesh data structure
#ifndef SWIG
std::vector<std::shared_ptr<Body>> bodies; ///< bodies in fluid
#endif
// Reference values
double alpha; ///< Angle of attack
double beta; ///< Angle of sideslip
+ double mach; ///< Mach number
double sR; ///< Reference surface
double cR; ///< Reference chord
Eigen::Vector3d xR; ///< Reference center point (for moment computation)
- Problem(std::shared_ptr<tbox::MshData> _smsh, double aoa, double aos, double sref, double cref, double xref, double yref, double zref);
+ Problem(std::shared_ptr<tbox::MshData> _msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref);
~Problem() { std::cout << "~Problem()\n"; }
void add(std::shared_ptr<Body> b);
+ // Functions
+ inline double PhiI(Eigen::Vector3d const &pos);
+ inline Eigen::Vector3d Ui();
+ Eigen::Vector3d U(tbox::Element &e, std::vector<double> const &phi);
+ double Rho(tbox::Element &e, std::vector<double> const &phi);
+ double Mach(tbox::Element &e, std::vector<double> const &phi);
+ double Cp(tbox::Element &e, std::vector<double> const &phi);
+
#ifndef SWIG
void check() const;
void updateMem();
@@ -55,6 +65,8 @@ public:
#endif
};
+#include "fProblem.inl"
+
} // namespace fpm
#endif //FPROBLEM_H
diff --git a/fpm/src/fProblem.inl b/fpm/src/fProblem.inl
new file mode 100644
index 0000000000000000000000000000000000000000..175e1a3d89ec95a0def6b9fba88769ec0cbe2f2c
--- /dev/null
+++ b/fpm/src/fProblem.inl
@@ -0,0 +1,31 @@
+/*
+ * 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.
+ */
+
+/**
+ * @brief Compute freestream potential
+ */
+inline double Problem::PhiI(Eigen::Vector3d const &pos)
+{
+ return cos(alpha) * cos(beta) * pos(0) + sin(beta) * pos(1) + sin(alpha) * cos(beta) * pos(2);
+}
+
+/**
+ * @brief Compute freestream velocity
+ */
+inline Eigen::Vector3d Problem::Ui()
+{
+ return Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta));
+}
diff --git a/fpm/src/fSolver.cpp b/fpm/src/fSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1769580d86d10db06ec63b1ee127bfe1b407fb97
--- /dev/null
+++ b/fpm/src/fSolver.cpp
@@ -0,0 +1,228 @@
+/*
+ * 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 "fSolver.h"
+#include "fBuilder.h"
+#include "fProblem.h"
+#include "fBody.h"
+
+#include "wMshData.h"
+#include "wNode.h"
+#include "wElement.h"
+#include "wMem.h"
+#include "wTag.h"
+#include "wResults.h"
+#include "wMshExport.h"
+#include "wCache.h"
+
+using namespace tbox;
+using namespace fpm;
+
+#define ANSI_COLOR_YELLOW "\x1b[1;33m"
+#define ANSI_COLOR_RESET "\x1b[0m"
+
+/**
+ * @brief Initialize the solver and perform sanity checks
+ * @authors Adrien Crovato
+ */
+Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0)
+{
+ // Say hi
+ std::cout << std::endl;
+ std::cout << "*******************************************************************************" << std::endl;
+ std::cout << "** \\_/ **" << std::endl;
+ std::cout << "** \\_______O(_)O_______/ **" << std::endl;
+ std::cout << "** **" << std::endl;
+ std::cout << "*******************************************************************************" << std::endl;
+ std::cout << "** Hi! My name is FPM v0.1-2007 **" << std::endl;
+ std::cout << "** Adrien Crovato & Romain Boman **" << std::endl;
+ std::cout << "** ULiege 2020-2021 **" << std::endl;
+ std::cout << "*******************************************************************************" << std::endl
+ << std::endl;
+
+ // Setup variables (unkwown vector)
+ phi.resize(aic->pbl->msh->nodes.size(), 0.);
+ vPhi.resize(aic->pbl->msh->nodes.size(), 0.);
+ U.resize(aic->pbl->msh->elems.size(), Eigen::Vector3d::Zero());
+ rho.resize(aic->pbl->msh->elems.size(), 0.);
+ M.resize(aic->pbl->msh->elems.size(), 0.);
+ Cp.resize(aic->pbl->msh->elems.size(), 0.);
+}
+
+/**
+ * @brief Solve the (full) potential equation
+ */
+void Solver::run()
+{
+ if (!aic->valid)
+ throw std::runtime_error("fpm::Solver:run() AIC matrices are not up-to-date!\n");
+
+ // @todo compute sigma from an interpolator and compute volume flow
+ std::cout << "Computing field sources... " << std::flush;
+ Eigen::VectorXd sigma = Eigen::VectorXd::Zero(aic->nF);
+ Eigen::VectorXd uX = aic->Cx * sigma;
+ Eigen::VectorXd uY = aic->Cy * sigma;
+ Eigen::VectorXd uZ = aic->Cz * sigma;
+ std::cout << "done!" << std::endl;
+
+ std::cout << "Computing body sources... " << std::flush;
+ Eigen::VectorXd mu(aic->nP), tau(aic->nP);
+ int ig = 0;
+ for (auto body : aic->pbl->bodies)
+ for (int i = 0; i < body->tag->elems.size(); ++i, ++ig)
+ tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal());
+ std::cout << "done!" << std::endl;
+
+ std::cout << "Solving for body doublets... " << std::flush;
+ Eigen::PartialPivLU<Eigen::MatrixXd>(aic->A).solve(aic->B * tau);
+ std::cout << "done!" << std::endl;
+
+ std::cout << "Computing flow and loads on bodies... " << std::flush;
+ computeFlow(mu, tau, sigma);
+ computeLoad();
+ std::cout << "done!" << std::endl << std::endl;
+}
+
+/**
+ * @brief Write the results
+ */
+void Solver::save(std::shared_ptr<MshExport> mshWriter, int n)
+{
+ // Write files
+ std::cout << "Saving files... " << std::endl;
+ // setup results
+ Results results;
+ results.scalars_at_nodes["phi"] = φ
+ results.scalars_at_nodes["varPhi"] = &vPhi;
+ results.vectors_at_elems["U"] = &U;
+ results.scalars_at_elems["rho"] = ρ
+ results.scalars_at_elems["Mach"] = &M;
+ results.scalars_at_elems["Cp"] = &Cp;
+ // save (mesh and boundary surface)
+ if (n > 0)
+ {
+ mshWriter->save(aic->pbl->msh->name + "_" + std::to_string(n), results);
+ for (auto body : aic->pbl->bodies)
+ body->save(body->tag->name + "_" + std::to_string(n));
+ }
+ else
+ {
+ mshWriter->save(aic->pbl->msh->name, results);
+ for (auto body : aic->pbl->bodies)
+ body->save(body->tag->name);
+ }
+}
+
+/**
+ * @brief Compute flow solution on bodies
+ */
+void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma)
+{
+ // @todo should we
+ // - compute velocity from phi at nodes or mu at nodes? if mu, then store phi at elements only
+
+ // Compute phi at elements and transfer at nodes
+ Eigen::VectorXd phie = aic->A * mu + aic->B * tau + aic->C * sigma;
+ auto pbl = aic->pbl;
+ // reset all values
+ for (auto n : pbl->msh->nodes)
+ vPhi[n->no - 1] = 0;
+ int ig = 0;
+ for (auto body : pbl->bodies)
+ {
+ // associate the potential value to its element
+ std::map<Element *, double> mapPhi;
+ for (auto e : body->tag->elems)
+ {
+ mapPhi[e] = phie[ig];
+ ig++;
+ }
+ // average at nodes
+ for (auto n : body->nodes)
+ for (auto e : body->neMap.at(n))
+ vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); // maybe use area average?
+ }
+
+ // Compute surface flow
+ for (auto body : pbl->bodies)
+ {
+ for (auto n : body->nodes)
+ phi[n->no - 1] = pbl->PhiI(n->pos) + vPhi[n->row]; // total potential
+
+ for (auto e : body->tag->elems)
+ {
+ U[e->no - 1] = pbl->U(*e, phi); // velocity
+ rho[e->no - 1] = pbl->Rho(*e, phi); // density
+ M[e->no - 1] = pbl->Mach(*e, phi); // Mach number
+ Cp[e->no - 1] = pbl->Cp(*e, phi); // pressure coefficient
+ }
+ }
+}
+
+/**
+ * @brief Compute aerodynamic load coefficients on bodies
+ */
+void Solver::computeLoad()
+{
+ Cl = 0;
+ Cd = 0;
+ Cs = 0;
+ Cm = 0;
+ auto pbl = aic->pbl;
+ for (auto body : pbl->bodies)
+ {
+ // Compute normalized load (i.e. load / pressure)
+ for (size_t i = 0; i < body->nodes.size(); ++i)
+ {
+ body->cLoadX[i] = 0;
+ body->cLoadY[i] = 0;
+ body->cLoadZ[i] = 0;
+ for (auto e : body->neMap.at(body->nodes[i]))
+ {
+ Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVMem().getVol() * e->normal() / e->nodes.size();
+ body->cLoadX[i] += cLoad(0);
+ body->cLoadY[i] += cLoad(1);
+ body->cLoadZ[i] += cLoad(2);
+ }
+ }
+
+ // Compute aerodynamic load coefficients (i.e. Load / pressure / surface)
+ // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
+ body->Cm = 0;
+ Eigen::Vector3d Cf(0, 0, 0);
+ for (auto e : body->tag->elems)
+ {
+ Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR; // lever arm
+ Eigen::Vector3d cf = Cp[e->no - 1] * e->getVMem().getVol() * e->normal(); // force coefficient
+ Cf -= cf;
+ body->Cm -= (l(2) * cf(0) - l(0) * cf(2)) / (pbl->sR * pbl->cR); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+ }
+ // rotate to flow direction
+ body->Cd = (Cf(0) * cos(pbl->alpha) * cos(pbl->beta) + Cf(1) * sin(pbl->beta) + Cf(2) * sin(pbl->alpha) * cos(pbl->beta)) / pbl->sR;
+ body->Cs = (-Cf(0) * cos(pbl->alpha) * sin(pbl->beta) + Cf(1) * cos(pbl->beta) - Cf(2) * sin(pbl->alpha) * sin(pbl->beta)) / pbl->sR;
+ body->Cl = (-Cf(0) * sin(pbl->alpha) + Cf(2) * cos(pbl->alpha)) / pbl->sR;
+ // compute total
+ Cl += body->Cl;
+ Cd += body->Cd;
+ Cs += body->Cs;
+ Cm += body->Cm;
+ }
+}
+
+void Solver::write(std::ostream &out) const
+{
+ out << "fpm::Solver\nwith\t" << *aic;
+}
diff --git a/fpm/src/fSolver.h b/fpm/src/fSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..a7c90841450d68adcddef23a7e8b6d15267fbe0e
--- /dev/null
+++ b/fpm/src/fSolver.h
@@ -0,0 +1,69 @@
+/*
+ * 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.
+ */
+
+#ifndef FSOLVER_H
+#define FSOLVER_H
+
+#include "fpm.h"
+#include "wObject.h"
+
+#include <iostream>
+#include <vector>
+#include <memory>
+#include <Eigen/Dense>
+
+namespace fpm
+{
+
+/**
+ * @brief Field panel solver
+ * @authors Adrien Crovato
+ * @todo check and implement velocity computation
+ */
+class FPM_API Solver : public fwk::wSharedObject
+{
+public:
+ std::shared_ptr<Builder> aic; ///< AIC builder
+
+ std::vector<double> phi; ///< full potential
+ std::vector<double> vPhi; ///< perturbation potential
+ std::vector<Eigen::Vector3d> U; ///< velocity
+ std::vector<double> rho; ///< density
+ std::vector<double> M; ///< mach number
+ std::vector<double> Cp; ///< pressure coefficient
+ double Cl; ///< lift coefficient
+ double Cd; ///< drag coefficient
+ double Cs; ///< sideforce coefficient
+ double Cm; ///< pitch moment coefficient (positive nose-up)
+
+public:
+ Solver(std::shared_ptr<Builder> _aic);
+ ~Solver() { std::cout << "~Solver()\n"; }
+
+ void run();
+ void save(std::shared_ptr<tbox::MshExport> mshWriter, int n = 0);
+
+#ifndef SWIG
+ virtual void write(std::ostream &out) const override;
+#endif
+
+private:
+ void computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma);
+ void computeLoad();
+};
+
+} // namespace fpm
+#endif //FSOLVER_H
diff --git a/fpm/src/fpm.h b/fpm/src/fpm.h
index a56c56666b8cef734fd365b5fdff0158dab358e1..b9b586efeb139c9cdd1fe7828948ea7e4cc2875e 100644
--- a/fpm/src/fpm.h
+++ b/fpm/src/fpm.h
@@ -43,6 +43,10 @@ class Problem;
class Body;
class Wake;
+// solver
+class Builder;
+class Solver;
+
// misc
class Edge;
class EquEdge;
diff --git a/fpm/tests/basic.py b/fpm/tests/basic.py
index af8a994832dc1d3edcabf841750934bbd74664c3..8228946bcb6a0a529840386a0aeacd235a2f3a61 100644
--- a/fpm/tests/basic.py
+++ b/fpm/tests/basic.py
@@ -33,10 +33,17 @@ def main():
pars = {'spn' : 5, 'nC' : 80, 'nS' : 10}
msh = gmsh.MeshLoader('../models/n0012.geo', __file__).execute(**pars)
# problem
- pbl = fpm.Problem(msh, 0, 0, 5., 1., 0., 0., 0.)
+ pbl = fpm.Problem(msh, 0, 0, 0, 5., 1., 0., 0., 0.)
wing = fpm.Body(msh, 'wing', ['te'], 5)
pbl.add(wing)
- print(pbl)
+ # AIC builder
+ aic = fpm.Builder(pbl)
+ aic.run()
+ # solver
+ slv = fpm.Solver(aic)
+ slv.run()
+ slv.save(tbox.GmshExport(msh))
+ print('\n', slv)
# end timer
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)