From 5ddb2c5a8d088bdfbfbd381c0a1d63273c000e48 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Mon, 14 Jun 2021 18:14:23 +0200
Subject: [PATCH] Add formulation classes to build set of equations, and remove
build methods from Element (flow only). Draft methods to compute freestream
velocity. Rework data structure of boundary. Add const. Improve test cases.
---
flow/src/flow.h | 8 +
flow/src/wAdjoint.cpp | 331 +++++++++----------------------
flow/src/wBlowingResidual.cpp | 48 +++++
flow/src/wBlowingResidual.h | 39 ++++
flow/src/wBoundary.cpp | 25 +--
flow/src/wBoundary.h | 21 +-
flow/src/wF0El.cpp | 24 +--
flow/src/wF0El.h | 24 +--
flow/src/wF1El.cpp | 15 +-
flow/src/wF1El.h | 10 +-
flow/src/wFreestreamResidual.cpp | 70 +++++++
flow/src/wFreestreamResidual.h | 41 ++++
flow/src/wKuttaElement.cpp | 164 ---------------
flow/src/wKuttaElement.h | 15 +-
flow/src/wKuttaResidual.cpp | 167 ++++++++++++++++
flow/src/wKuttaResidual.h | 44 ++++
flow/src/wLoadFunctional.cpp | 136 +++++++++++++
flow/src/wLoadFunctional.h | 42 ++++
flow/src/wNewton.cpp | 103 +++-------
flow/src/wPicard.cpp | 19 +-
flow/src/wPotentialResidual.cpp | 250 +++++++++++++++++++++++
flow/src/wPotentialResidual.h | 45 +++++
flow/src/wProblem.cpp | 62 ++++++
flow/src/wProblem.h | 3 +
flow/src/wSolver.cpp | 59 +++---
flow/src/wWakeElement.cpp | 293 +--------------------------
flow/src/wWakeElement.h | 17 +-
flow/src/wWakeResidual.cpp | 251 +++++++++++++++++++++++
flow/src/wWakeResidual.h | 44 ++++
flow/tests/lift.py | 3 +
flow/tests/lift3.py | 2 +
flow/tests/meshDef.py | 4 +-
flow/tests/meshDef3.py | 2 +
flow/tests/nonlift.py | 3 +-
tbox/src/wElement.cpp | 12 +-
tbox/src/wElement.h | 6 +-
tbox/src/wFct0.cpp | 6 +-
tbox/src/wFct0.h | 12 +-
tbox/src/wFct1.h | 13 +-
tbox/src/wFct2.cpp | 8 +-
tbox/src/wFct2.h | 12 +-
tbox/src/wLine2.cpp | 4 +-
tbox/src/wLine2.h | 2 +-
tbox/src/wMem.inl | 2 +-
tbox/src/wTetra4.cpp | 48 +----
tbox/src/wTetra4.h | 6 +-
tbox/src/wTri3.cpp | 48 +----
tbox/src/wTri3.h | 6 +-
48 files changed, 1542 insertions(+), 1027 deletions(-)
create mode 100644 flow/src/wBlowingResidual.cpp
create mode 100644 flow/src/wBlowingResidual.h
create mode 100644 flow/src/wFreestreamResidual.cpp
create mode 100644 flow/src/wFreestreamResidual.h
create mode 100644 flow/src/wKuttaResidual.cpp
create mode 100644 flow/src/wKuttaResidual.h
create mode 100644 flow/src/wLoadFunctional.cpp
create mode 100644 flow/src/wLoadFunctional.h
create mode 100644 flow/src/wPotentialResidual.cpp
create mode 100644 flow/src/wPotentialResidual.h
create mode 100644 flow/src/wWakeResidual.cpp
create mode 100644 flow/src/wWakeResidual.h
diff --git a/flow/src/flow.h b/flow/src/flow.h
index 20e132f0..351ca549 100644
--- a/flow/src/flow.h
+++ b/flow/src/flow.h
@@ -51,6 +51,14 @@ class KuttaElement;
class Blowing;
class Boundary;
+// Formulation
+class PotentialResidual;
+class FreestreamResidual;
+class WakeResidual;
+class KuttaResidual;
+class BlowingResidual;
+class LoadFunctional;
+
// Solver
class Solver;
class Picard;
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index 225327f0..6774d1bb 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -22,12 +22,15 @@
#include "wKutta.h"
#include "wBoundary.h"
#include "wNewton.h"
+#include "wPotentialResidual.h"
+#include "wFreestreamResidual.h"
+#include "wWakeResidual.h"
+#include "wKuttaResidual.h"
+#include "wLoadFunctional.h"
#include "wMshData.h"
#include "wNode.h"
#include "wElement.h"
-#include "wCache.h"
-#include "wMem.h"
#include "wLinearSolver.h"
#include "wMshDeform.h"
#include "wResults.h"
@@ -74,11 +77,10 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
}
/**
- * @brief Run the linear solver
+ * @brief Run the adjoint method
*
* Solve the adjoint steady transonic Full Potential Equation
- * i.e. computes sensitivites for lift and drag cost functions
- * @authors Adrien Crovato
+ * and compute lift and drag sensitivites
*/
void Adjoint::run()
{
@@ -149,6 +151,7 @@ void Adjoint::solve()
/**
* @brief Compute total gradients of the loads with respect to angle of attack
+ * @todo one gradient for each surface node
*/
void Adjoint::computeGradientAoa()
{
@@ -237,6 +240,7 @@ void Adjoint::computeGradientMesh()
/**
* @brief Build the gradients of the lift and drag with respect to the flow variables
+ * @todo one gradient for each surface node
*/
void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
{
@@ -247,32 +251,30 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
dL.setZero();
dD.setZero();
+ // Lift and drag normalized directions
+ Eigen::RowVector3d dirL, dirD;
+ std::tie(dirD, std::ignore, dirL) = sol->pbl->computeDirections();
+
// Gradient of lift and drag
for (auto sur : sol->pbl->bnds)
{
tbb::parallel_do(sur->groups[0]->tag->elems.begin(), sur->groups[0]->tag->elems.end(), [&](Element *e) {
- // Build flux vector b = dCp * V * [inv(J)*dN]
+ // Associated volume element
Element *eV = sur->svMap.at(e);
- Eigen::RowVectorXd be = sol->pbl->medium->cP.evalD(*eV, sol->phi, 0) * eV->computeGrad(sol->phi, 0).transpose() * eV->getVMem().getJinv(0) * eV->getVCache().getDsf(0) * e->getVMem().getVol();
- // Rotate to flow direction
- double dirL, dirD;
- if (sol->pbl->nDim == 2)
- {
- dirL = -e->normal().dot(Eigen::Vector3d(-sin(sol->pbl->alpha), cos(sol->pbl->alpha), 0)) / sol->pbl->S_ref;
- dirD = -e->normal().dot(Eigen::Vector3d(cos(sol->pbl->alpha), sin(sol->pbl->alpha), 0)) / sol->pbl->S_ref;
- }
- else
- {
- dirL = -e->normal().dot(Eigen::Vector3d(-sin(sol->pbl->alpha), 0, cos(sol->pbl->alpha))) / sol->pbl->S_ref;
- dirD = -e->normal().dot(Eigen::Vector3d(cos(sol->pbl->alpha) * cos(sol->pbl->beta), sin(sol->pbl->beta), sin(sol->pbl->alpha) * cos(sol->pbl->beta))) / sol->pbl->S_ref;
- }
+ // Build
+ Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, sol->pbl->medium->cP);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < eV->nodes.size(); ++ii)
+ for (size_t i = 0; i < e->nodes.size(); ++i)
{
- Node *nodi = eV->nodes[ii];
- dL(nodi->row) += dirL * be(ii);
- dD(nodi->row) += dirD * be(ii);
+ Eigen::RowVectorXd AeL = -dirL * Ae.block(i * 3, 0, 3, eV->nodes.size());
+ Eigen::RowVectorXd AeD = -dirD * Ae.block(i * 3, 0, 3, eV->nodes.size());
+ for (size_t j = 0; j < eV->nodes.size(); ++j)
+ {
+ Node *nodj = eV->nodes[j];
+ dL(nodj->row) += AeL(j);
+ dD(nodj->row) += AeD(j);
+ }
}
});
}
@@ -289,24 +291,12 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
// Reset gradient
dR.setZero();
- // Gradient of freestream velocity
- Eigen::Vector3d dvi;
- if (sol->pbl->nDim == 2)
- dvi = Eigen::Vector3d(-sin(sol->pbl->alpha), cos(sol->pbl->alpha), 0);
- else
- dvi = Eigen::Vector3d(-sin(sol->pbl->alpha) * cos(sol->pbl->beta), 0, cos(sol->pbl->alpha) * cos(sol->pbl->beta));
-
// Freestream contribution
for (auto nBC : sol->pbl->fBCs)
{
tbb::parallel_do(nBC->tag->elems.begin(), nBC->tag->elems.end(), [&](Element *e) {
// Build elementary flux vector
- Cache &cache = e->getVCache();
- Gauss &gauss = cache.getVGauss();
- Mem &mem = e->getVMem();
- Eigen::VectorXd be = Eigen::VectorXd::Zero(e->nodes.size());
- for (size_t k = 0; k < gauss.getN(); ++k)
- be += cache.getSf(k) * dvi.dot(e->normal()) * gauss.getW(k) * mem.getDetJ(k);
+ Eigen::VectorXd be = FreestreamResidual::buildGradientAoa(*e, sol->phi, *nBC);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -322,10 +312,11 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
// @todo very slight influence, can be removed by using Vi = [1,0,0]. Same results, same cost, less general but less code
for (auto wake : sol->pbl->wBCs)
{
+ Eigen::VectorXd dvi = sol->pbl->computeGradientVi(); // gradient of freestream velocity
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build bu, bl
Eigen::VectorXd bue, ble;
- we->buildGradientResidualAoA(bue, ble, sol->phi, dvi.block(0, 0, sol->pbl->nDim, 1));
+ std::tie(bue, ble) = WakeResidual::build(*we, sol->phi, dvi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -340,33 +331,25 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
/**
* @brief Compute gradients of the lift and drag with respect to angle of attack
- * @todo compute one gradient per boundary
+ * @todo one gradient for each surface node
*/
void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
{
- // Reste gradients
- dL = 0, dD = 0;
- // Compute aerodynamic load coefficients (i.e. Load / pressure / surface) along x and vertical (y in 2D, z in 3D) directions
+ // Gradient of lift and drag normalized directions
+ Eigen::RowVector3d dirL, dirD;
+ std::tie(dirD, std::ignore, dirL) = sol->pbl->computeGradientDirections();
+ // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
Eigen::Vector3d Cf(0, 0, 0);
for (auto bnd : sol->pbl->bnds)
- for (auto e : bnd->groups[0]->tag->elems)
- Cf -= bnd->integrate(e, sol->phi, sol->pbl->medium->cP) * e->normal();
+ for (size_t i = 0; i < bnd->nodes.size(); ++i)
+ Cf -= bnd->nLoads[i];
// Rotate to gradient of flow direction
- if (sol->pbl->nDim == 2)
- {
- dD += (-Cf(0) * sin(sol->pbl->alpha) + Cf(1) * cos(sol->pbl->alpha)) / sol->pbl->S_ref;
- dL += (-Cf(0) * cos(sol->pbl->alpha) - Cf(1) * sin(sol->pbl->alpha)) / sol->pbl->S_ref;
- }
- else
- {
- dD += (-Cf(0) * sin(sol->pbl->alpha) * cos(sol->pbl->beta) + Cf(2) * cos(sol->pbl->alpha) * cos(sol->pbl->beta)) / sol->pbl->S_ref;
- dL += (-Cf(0) * cos(sol->pbl->alpha) - Cf(2) * sin(sol->pbl->alpha)) / sol->pbl->S_ref;
- }
+ dD = Cf.dot(dirD);
+ dL = Cf.dot(dirL);
}
/**
* @brief Compute gradients of the residual with respect to mesh nodes
- * @todo compute constant variables outside GP loop
*/
void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
{
@@ -383,130 +366,31 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
Element *e = p.first;
//Upwind element
Element *eU = p.second[0];
+ // Build elementary matrices
+ Eigen::MatrixXd Ae1, Ae2;
+ std::tie(Ae1, Ae2) = PotentialResidual::buildGradientMesh(*e, *eU, sol->phi, *fluid, sol->pbl->nDim, sol->muCv, sol->mCOv);
- // Get Jacobian, shape functions and gauss points
- Cache &cache = e->getVCache();
- Mem &mem = e->getVMem();
- Gauss &gauss = cache.getVGauss();
-
- // Subsonic contribution
- Eigen::MatrixXd Ae1 = Eigen::MatrixXd::Zero(e->nodes.size(), sol->pbl->nDim * e->nodes.size());
- for (size_t k = 0; k < gauss.getN(); ++k)
+ // Assembly (subsonic)
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t i = 0; i < e->nodes.size(); ++i)
{
- // @todo, use first GP only? -> better: have 1GP per element?
- Eigen::VectorXd dPhi = e->computeGrad(sol->phi, k);
- double rho = fluid->rho.eval(*e, sol->phi, k);
- double dRho = fluid->rho.evalD(*e, sol->phi, k);
- Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- // Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- std::vector<double> const &dDetJ = mem.getGradDetJ(k);
- //
- double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
- for (size_t i = 0; i < e->nodes.size(); ++i)
+ size_t rowi = e->nodes[i]->row;
+ for (size_t j = 0; j < e->nodes.size(); ++j)
{
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
+ for (int n = 0; n < sol->pbl->nDim; ++n)
{
- size_t idx = i * sol->pbl->nDim + m;
- Ae1.col(idx) += w * dRho * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // drho * grad_phi*grad_psi * detj
- Ae1.col(idx) += w * rho * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // rho * dgrad_phi * grad_psi * detj
- Ae1.col(idx) += w * rho * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // rho * grad_phi * dgrad_psi * detj
- Ae1.col(idx) += w * rho * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // rho * grad_phi * grad_psi * ddetj
+ for (size_t jj = 0; jj < arows[e->nodes[j]->row].size(); ++jj)
+ {
+ size_t rowj = sol->pbl->nDim * (arows[e->nodes[j]->row][jj]) + n;
+ T.push_back(Eigen::Triplet<double>(sol->rows[rowi], rowj, Ae1(i, sol->pbl->nDim * j + n)));
+ }
}
}
}
- // Supersonic contribution
+ // Assembly (supersonic)
double mach = fluid->mach.eval(*e, sol->phi, 0);
if (mach > sol->mCOv)
{
- // switching function and gradient
- double mu = sol->muCv * (1 - (sol->mCOv * sol->mCOv) / (mach * mach));
- double dmu = 2 * sol->muCv * sol->mCOv * sol->mCOv / (mach * mach * mach);
- // upwind density and gradient
- Eigen::VectorXd dPhiU = eU->computeGrad(sol->phi, 0);
- double rhoU = fluid->rho.eval(*eU, sol->phi, 0);
- double dRhoU = fluid->rho.evalD(*eU, sol->phi, 0);
- Eigen::MatrixXd const &iJU = eU->getVMem().getJinv(0);
-
- // scale subsonic terms
- Ae1 *= 1 - mu;
-
- // compute upwind part of subsonic terms
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // @todo, use first GP only? -> better: have 1GP per element?
- Eigen::VectorXd dPhi = e->computeGrad(sol->phi, k);
- Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- // Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- std::vector<double> const &dDetJ = mem.getGradDetJ(k);
- //
- double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
- for (size_t i = 0; i < e->nodes.size(); ++i)
- {
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
- {
- size_t idx = i * sol->pbl->nDim + m;
- Ae1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // mu * rhoU * dgrad_phi * grad_psi * detj
- Ae1.col(idx) += w * mu * rhoU * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // mu * rhoU * grad_phi * dgrad_psi * detj
- Ae1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // mu * rhoU * grad_phi * grad_psi * ddetj
- }
- }
- }
-
- // compute derivative of switching function terms
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // @todo, use first GP only? -> better: have 1GP per element!
- Eigen::VectorXd dPhi = e->computeGrad(sol->phi, k);
- double rho = fluid->rho.eval(*e, sol->phi, k);
- double dM = fluid->mach.evalD(*e, sol->phi, k);
- Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- // Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- //
- double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
- for (size_t i = 0; i < e->nodes.size(); ++i)
- {
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
- {
- size_t idx = i * sol->pbl->nDim + m;
- Ae1.col(idx) += w * (-rho + rhoU) * dmu * dM * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // (rhoU-rho) * dmu * grad_phi * grad_psi * detj
- }
- }
- }
-
- // compute derivative of upwind density terms (on upwind element)
- Eigen::MatrixXd Ae2 = Eigen::MatrixXd::Zero(e->nodes.size(), sol->pbl->nDim * eU->nodes.size());
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // @todo, use first GP only? -> better: have 1GP per element!
- Eigen::VectorXd dPhi = e->computeGrad(sol->phi, k);
- Eigen::MatrixXd const &dSf = cache.getDsf(k);
- Eigen::MatrixXd const &iJ = mem.getJinv(k);
- // Jacobian gradients
- std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
- //
- double w = gauss.getW(k);
- double detJ = mem.getDetJ(k);
- for (size_t i = 0; i < eU->nodes.size(); ++i)
- {
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
- {
- size_t idx = i * sol->pbl->nDim + m;
- Ae2.col(idx) += w * mu * dRhoU * dPhiU.transpose() * (-iJU * dJ[idx] * dPhiU) * (iJ * dSf).transpose() * dPhi * detJ; // mu * drhoU * grad_phi * grad_psi * detj
- }
- }
- }
-
- // Assembly (supersonic)
- tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
size_t rowi = e->nodes[i]->row;
@@ -523,23 +407,6 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
}
}
}
- // Assembly (subsonic)
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t i = 0; i < e->nodes.size(); ++i)
- {
- size_t rowi = e->nodes[i]->row;
- for (size_t j = 0; j < e->nodes.size(); ++j)
- {
- for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[e->nodes[j]->row].size(); ++jj)
- {
- size_t rowj = sol->pbl->nDim * (arows[e->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(sol->rows[rowi], rowj, Ae1(i, sol->pbl->nDim * j + n)));
- }
- }
- }
- }
});
// Farfield B.C. are not taken into account because they have no influence on the solution
@@ -550,7 +417,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build elementary matrices
Eigen::MatrixXd Aeu, Ael, Aes;
- we->buildGradientResidualMesh(Aeu, Ael, Aes, sol->phi, vi);
+ std::tie(Aeu, Ael, Aes) = WakeResidual::buildGradientMesh(*we, sol->phi, vi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < we->nRow; ++i)
@@ -598,7 +465,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
tbb::parallel_do(kutta->kEle.begin(), kutta->kEle.end(), [&](KuttaElement *ke) {
// Build elementary matrices
Eigen::MatrixXd Ae, Aes;
- ke->buildGradientResidualMesh(Ae, Aes, sol->phi);
+ std::tie(Ae, Aes) = KuttaResidual::buildGradientMesh(*ke, sol->phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < ke->nRow; ++i)
@@ -637,7 +504,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
/**
* @brief Compute gradients of the lift and drag with respect to mesh nodes
- * @todo compute one gradient per boundary
+ * @todo one gradient for each surface node
*/
void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd &dD)
{
@@ -648,68 +515,54 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd
dL.setZero();
dD.setZero();
- Eigen::MatrixXd dCf = Eigen::MatrixXd::Zero(3, dL.size());
+ // Lift and drag normalized directions
+ Eigen::RowVector3d dirL, dirD;
+ std::tie(dirD, std::ignore, dirL) = sol->pbl->computeDirections();
+
+ // Gradient of lift and drag
for (auto bnd : sol->pbl->bnds)
{
tbb::parallel_do(bnd->groups[0]->tag->elems.begin(), bnd->groups[0]->tag->elems.end(), [&](Element *e) {
// Volume element
Element *eV = bnd->svMap.at(e);
- // Variables
- double cp = sol->pbl->medium->cP.eval(*eV, sol->phi, 0); // pressure coefficient
- Eigen::Vector3d nrm = e->normal(); // unit normal
- double surf = e->getVMem().getVol(); // surface
- // Gradients
- double dCp = sol->pbl->medium->cP.evalD(*eV, sol->phi, 0); // pressure coefficient
- Eigen::VectorXd dPhi = eV->computeGrad(sol->phi, 0); // potential
- Eigen::MatrixXd const &iJ = eV->getVMem().getJinv(0); // inverse Jacobian
- std::vector<Eigen::MatrixXd> const &dJ = eV->getVMem().getGradJ(0); // Jacobian
- std::vector<Eigen::Vector3d> dNrm = e->gradientNormal(); // unit normal
- std::vector<double> dSurf = e->getVMem().getGradVol(); // surface
-
- // Build element contributions
- Eigen::MatrixXd dCfV = Eigen::MatrixXd::Zero(3, eV->nodes.size() * sol->pbl->nDim);
- Eigen::MatrixXd dCfS = Eigen::MatrixXd::Zero(3, e->nodes.size() * sol->pbl->nDim);
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
- {
- // volume
- for (size_t i = 0; i < eV->nodes.size(); ++i)
- {
- size_t idx = i * sol->pbl->nDim + m;
- dCfV.col(idx) = dCp * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * nrm * surf; // dCp * n * S
- }
- // surface
- for (size_t i = 0; i < e->nodes.size(); ++i)
- {
- size_t idx = i * sol->pbl->nDim + m;
- dCfS.col(idx) = cp * dNrm[idx] * surf + cp * nrm * dSurf[idx]; // Cp * dn * S + Cp * n * dS
- }
- }
+ // Build
+ Eigen::MatrixXd Aes, Aev;
+ std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, sol->pbl->medium->cP, sol->pbl->nDim);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t m = 0; m < sol->pbl->nDim; ++m)
+ for (size_t i = 0; i < e->nodes.size(); ++i)
{
- // volume
- for (size_t i = 0; i < eV->nodes.size(); ++i)
- for (size_t ii = 0; ii < arows[eV->nodes[i]->row].size(); ++ii)
- dCf.col((arows[eV->nodes[i]->row][ii]) * sol->pbl->nDim + m) -= dCfV.col(i * sol->pbl->nDim + m);
- // surface
- for (size_t i = 0; i < e->nodes.size(); ++i)
- for (size_t ii = 0; ii < arows[e->nodes[i]->row].size(); ++ii)
- dCf.col((arows[e->nodes[i]->row][ii]) * sol->pbl->nDim + m) -= dCfS.col(i * sol->pbl->nDim + m);
+ // rotate to flow direction
+ Eigen::RowVectorXd AesL = -dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AesD = -dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AevL = -dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+ Eigen::RowVectorXd AevD = -dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+ for (size_t m = 0; m < sol->pbl->nDim; ++m)
+ {
+ // surface
+ for (size_t j = 0; j < e->nodes.size(); ++j)
+ {
+ size_t rowj = e->nodes[j]->row;
+ for (size_t jj = 0; jj < arows[rowj].size(); ++jj)
+ {
+ dL.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AesL.col(j * sol->pbl->nDim + m);
+ dD.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AesD.col(j * sol->pbl->nDim + m);
+ }
+ }
+ // volume
+ for (size_t j = 0; j < eV->nodes.size(); ++j)
+ {
+ size_t rowj = eV->nodes[j]->row;
+ for (size_t jj = 0; jj < arows[rowj].size(); ++jj)
+ {
+ dL.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AevL.col(j * sol->pbl->nDim + m);
+ dD.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AevD.col(j * sol->pbl->nDim + m);
+ }
+ }
+ }
}
});
}
- // Rotate to flow direction
- if (sol->pbl->nDim == 2)
- {
- dD = Eigen::RowVector3d(cos(sol->pbl->alpha), sin(sol->pbl->alpha), 0.) * dCf / sol->pbl->S_ref;
- dL = Eigen::RowVector3d(-sin(sol->pbl->alpha), cos(sol->pbl->alpha), 0.) * dCf / sol->pbl->S_ref;
- }
- else
- {
- dD = Eigen::RowVector3d(cos(sol->pbl->alpha) * cos(sol->pbl->beta), sin(sol->pbl->beta), sin(sol->pbl->alpha) * cos(sol->pbl->beta)) * dCf / sol->pbl->S_ref;
- dL = Eigen::RowVector3d(-sin(sol->pbl->alpha), 0., cos(sol->pbl->alpha)) * dCf / sol->pbl->S_ref;
- }
}
/**
diff --git a/flow/src/wBlowingResidual.cpp b/flow/src/wBlowingResidual.cpp
new file mode 100644
index 00000000..09dfcf30
--- /dev/null
+++ b/flow/src/wBlowingResidual.cpp
@@ -0,0 +1,48 @@
+/*
+ * 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 "wBlowingResidual.h"
+#include "wFct0.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Build the residual vector, on one boundary element
+ *
+ * b = \int psi * vn dV
+ */
+Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, Fct0 const &f)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Blowing velocity
+ double vn = f.eval(e, phi, 0);
+
+ // Integral of (normal) transpiration velocity
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+ return b;
+}
diff --git a/flow/src/wBlowingResidual.h b/flow/src/wBlowingResidual.h
new file mode 100644
index 00000000..95e09a3e
--- /dev/null
+++ b/flow/src/wBlowingResidual.h
@@ -0,0 +1,39 @@
+/*
+ * 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 WBLOWINGRESIDUAL_H
+#define WBLOWINGRESIDUAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of transpiration boundary terms
+ */
+class FLOW_API BlowingResidual
+{
+public:
+ // Newton
+ static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, tbox::Fct0 const &f);
+};
+
+} // namespace flow
+#endif //WBLOWINGRESIDUAL_H
diff --git a/flow/src/wBoundary.cpp b/flow/src/wBoundary.cpp
index ff6822b4..61e3bfb1 100644
--- a/flow/src/wBoundary.cpp
+++ b/flow/src/wBoundary.cpp
@@ -68,11 +68,9 @@ void Boundary::create()
std::sort(nodes.begin(), nodes.end());
auto iterator = std::unique(nodes.begin(), nodes.end());
nodes.resize(std::distance(nodes.begin(), iterator));
-
- // Create a map between the boundary (surface) nodes and elements
- for (auto e : groups[0]->tag->elems)
- for (auto n : e->nodes)
- neMap[n].push_back(e);
+ // Create a map between the boundary (surface) nodes and their local index
+ for (size_t i = 0; i < nodes.size(); ++i)
+ nMap[nodes[i]] = i;
// Create unordered set of surface faces
std::unordered_set<Face *, EquFace, EquFace> bndFaces;
@@ -113,22 +111,9 @@ void Boundary::create()
svMap[f->el0] = f->el1;
delete f;
}
- // Resize load vectors
- cLoadX.resize(nodes.size());
- cLoadY.resize(nodes.size());
- cLoadZ.resize(nodes.size());
-}
-/*
- * @brief Integrate f over an element
- * @authors Adrien Crovato
- */
-double Boundary::integrate(Element *e, std::vector<double> const &u, tbox::Fct0 const &f) const
-{
- // evaluate the function in the volume element
- double fk = f.eval(*svMap.at(e), u, 0);
- // integrate the function on the surface element
- return fk * e->getVMem().getVol(); // since fk is constant: sum(fk*wk*dtmk) = fk*V
+ // Resize loads vector
+ nLoads.resize(nodes.size(), Eigen::Vector3d::Zero());
}
/**
diff --git a/flow/src/wBoundary.h b/flow/src/wBoundary.h
index b40e70c1..1d4a1285 100644
--- a/flow/src/wBoundary.h
+++ b/flow/src/wBoundary.h
@@ -19,8 +19,10 @@
#include "flow.h"
#include "wGroups.h"
+
#include <vector>
#include <map>
+#include <Eigen/Dense>
namespace flow
{
@@ -32,22 +34,19 @@ namespace flow
class FLOW_API Boundary : public tbox::Groups
{
public:
- std::vector<tbox::Node *> nodes; ///< nodes of Boundary
- std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< node to element map
- std::map<tbox::Element *, tbox::Element *> svMap; ///< surface to volume element map
- 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<tbox::Node *> nodes; ///< nodes of Boundary
+ std::map<tbox::Node *, size_t> nMap; ///< node to local index map
+ std::map<tbox::Element *, tbox::Element *> svMap; ///< surface to volume element map
+ double Cl; ///< lift coefficient
+ double Cd; ///< drag coefficient
+ double Cs; ///< sideforce coefficient
+ double Cm; ///< pitch moment coefficient (positive nose-up)
+ std::vector<Eigen::Vector3d> nLoads; ///< nodal (normalized) aerodynamic loads on boundary
Boundary(std::shared_ptr<tbox::MshData> _msh, std::vector<int> const &nos);
Boundary(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
virtual ~Boundary() { std::cout << "~Boundary()\n"; }
- double integrate(tbox::Element *e, std::vector<double> const &u, tbox::Fct0 const &f) const;
void save(std::string const &name, tbox::Results const &res);
#ifndef SWIG
diff --git a/flow/src/wF0El.cpp b/flow/src/wF0El.cpp
index f4222fbd..45a27a4c 100644
--- a/flow/src/wF0El.cpp
+++ b/flow/src/wF0El.cpp
@@ -22,7 +22,7 @@ using namespace flow;
* @brief Evaluate the nonlinear density (constant over element)
* @authors Adrien Crovato
*/
-double F0ElRho::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU = e.computeGrad(u, k).norm(); // norm of velocity
@@ -43,7 +43,7 @@ double F0ElRho::eval(Element &e, std::vector<double> const &u, size_t k) const
* @brief Evaluate the nonlinear density derivative (constant over element)
* @authors Adrien Crovato
*/
-double F0ElRho::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU = e.computeGrad(u, k).norm(); // norm of velocity
@@ -69,7 +69,7 @@ void F0ElRho::write(std::ostream &out) const
* @brief Evaluate the linear density (constant over element)
* @authors Adrien Crovato
*/
-double F0ElRhoL::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
return 1.;
}
@@ -77,7 +77,7 @@ double F0ElRhoL::eval(Element &e, std::vector<double> const &u, size_t k) const
* @brief Evaluate the linear density derivative (constant over element)
* @authors Adrien Crovato
*/
-double F0ElRhoL::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
return 0.;
}
@@ -90,7 +90,7 @@ void F0ElRhoL::write(std::ostream &out) const
* @brief Evaluate the nonlinear Mach number (constant over element)
* @authors Adrien Crovato
*/
-double F0ElMach::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU2 = e.computeGrad2(u, k); // velocity squared
double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2); // speed of sound squared
@@ -101,7 +101,7 @@ double F0ElMach::eval(Element &e, std::vector<double> const &u, size_t k) const
* @brief Evaluate the nonlinear Mach number derivative (constant over element)
* @authors Adrien Crovato
*/
-double F0ElMach::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU2 = e.computeGrad2(u, k);
double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2);
@@ -117,7 +117,7 @@ void F0ElMach::write(std::ostream &out) const
* @brief Evaluate the linear Mach number (constant over element)
* @authors Adrien Crovato
*/
-double F0ElMachL::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
return 0.;
}
@@ -125,7 +125,7 @@ double F0ElMachL::eval(Element &e, std::vector<double> const &u, size_t k) const
* @brief Evaluate the linear Mach number derivative (constant over element)
* @authors Adrien Crovato
*/
-double F0ElMachL::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElMachL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
return 0.;
}
@@ -138,7 +138,7 @@ void F0ElMachL::write(std::ostream &out) const
* @brief Evaluate the nonlinear pressure coefficient (constant over element)
* @authors Adrien Crovato
*/
-double F0ElCp::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU2 = e.computeGrad2(u, k); // velocity squared
//particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1);
@@ -153,7 +153,7 @@ double F0ElCp::eval(Element &e, std::vector<double> const &u, size_t k) const
* @brief Evaluate the nonlinear pressure coefficient derivative (constant over element)
* @authors Adrien Crovato
*/
-double F0ElCp::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU2 = e.computeGrad2(u, k); // velocity squared
//particularized: -2 * pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
@@ -173,11 +173,11 @@ void F0ElCp::write(std::ostream &out) const
* @brief Evaluate the linear pressure coefficient (constant over element)
* @authors Adrien Crovato
*/
-double F0ElCpL::eval(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
return 1 - e.computeGrad2(u, k);
}
-double F0ElCpL::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double F0ElCpL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
return -2.; // has to be multiplied by grad u to recover true derivative
}
diff --git a/flow/src/wF0El.h b/flow/src/wF0El.h
index 02155449..a92e4d4b 100644
--- a/flow/src/wF0El.h
+++ b/flow/src/wF0El.h
@@ -46,10 +46,10 @@ public:
}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -66,10 +66,10 @@ public:
F0ElRhoL() : Fct0() {}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -91,10 +91,10 @@ public:
F0ElMach(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -111,10 +111,10 @@ public:
F0ElMachL() : Fct0() {}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -136,10 +136,10 @@ public:
F0ElCp(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -157,10 +157,10 @@ public:
F0ElCpL() : Fct0() {}
#ifndef SWIG
- virtual double eval(tbox::Element &e,
+ virtual double eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
- virtual double evalD(tbox::Element &e,
+ virtual double evalD(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
diff --git a/flow/src/wF1El.cpp b/flow/src/wF1El.cpp
index 6cba8d94..a65f35c8 100644
--- a/flow/src/wF1El.cpp
+++ b/flow/src/wF1El.cpp
@@ -39,18 +39,31 @@ F1ElVi::F1ElVi(int _nDim, double alpha, double beta) : Fct1(), nDim(_nDim)
void F1ElVi::update(double alpha, double beta)
{
if (nDim == 2)
+ {
v << cos(alpha), sin(alpha), 0;
+ dv << -sin(alpha), cos(alpha), 0;
+ }
else
+ {
v << cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta);
+ dv << -sin(alpha) * cos(beta), 0, cos(alpha) * cos(beta);
+ }
}
/**
* @brief Return the freestream velocity
* @authors Adrien Crovato
*/
-Eigen::Vector3d F1ElVi::eval(Element &e, std::vector<double> const &u, size_t k) const
+Eigen::Vector3d F1ElVi::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
return v;
}
+/**
+ * @brief Return the gradient of freestream velocity with respect to angle of attack
+ */
+Eigen::Vector3d F1ElVi::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+{
+ return dv;
+}
void F1ElVi::write(std::ostream &out) const
{
out << "F1ElVi";
diff --git a/flow/src/wF1El.h b/flow/src/wF1El.h
index 065dc22b..7d869712 100644
--- a/flow/src/wF1El.h
+++ b/flow/src/wF1El.h
@@ -31,17 +31,21 @@ namespace flow
*/
class FLOW_API F1ElVi : public tbox::Fct1
{
- int nDim; ///< number of dimensions
- Eigen::Vector3d v; ///< freestream velocity vector
+ int nDim; ///< number of dimensions
+ Eigen::Vector3d v; ///< freestream velocity vector
+ Eigen::Vector3d dv; ///< freestream velocity vector gradient (wrt alpha)
public:
F1ElVi(int _nDim, double alpha, double beta = 0);
void update(double alpha, double beta = 0);
#ifndef SWIG
- virtual Eigen::Vector3d eval(tbox::Element &e,
+ virtual Eigen::Vector3d eval(tbox::Element const &e,
std::vector<double> const &u,
size_t k) const override;
+ virtual Eigen::Vector3d evalD(tbox::Element const &e,
+ std::vector<double> const &u,
+ size_t k) const override;
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/flow/src/wFreestreamResidual.cpp b/flow/src/wFreestreamResidual.cpp
new file mode 100644
index 00000000..afdd43ae
--- /dev/null
+++ b/flow/src/wFreestreamResidual.cpp
@@ -0,0 +1,70 @@
+/*
+ * 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 "wFreestreamResidual.h"
+#include "wFreestream.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Build the residual vector, on one freestream boundary element
+ *
+ * b = \int psi * v_inf * n dV
+ */
+Eigen::VectorXd FreestreamResidual::build(Element const &e, std::vector<double> const &phi, Freestream const &fs)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Normal velocity
+ double vn = fs.f.eval(e, phi, 0).dot(e.normal());
+
+ // Integral of mass flux (density times normal velocity) through freestream boundary
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+ return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the angle of attack, on one freestream boundary element
+ *
+ * b = \int psi * dv_inf * n dV
+ */
+Eigen::VectorXd FreestreamResidual::buildGradientAoa(tbox::Element const &e, std::vector<double> const &phi, Freestream const &fs)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Gradient of normal velocity
+ double dvn = fs.f.evalD(e, phi, 0).dot(e.normal());
+
+ // Build mass flux gradient
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ b += cache.getSf(k) * dvn * gauss.getW(k) * mem.getDetJ(k);
+ return b;
+}
diff --git a/flow/src/wFreestreamResidual.h b/flow/src/wFreestreamResidual.h
new file mode 100644
index 00000000..abaa783e
--- /dev/null
+++ b/flow/src/wFreestreamResidual.h
@@ -0,0 +1,41 @@
+/*
+ * 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 WFREESTREAMRESIDUAL_H
+#define WFREESTREAMRESIDUAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of freestream boundary residuals
+ */
+class FLOW_API FreestreamResidual
+{
+public:
+ // Newton
+ static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, Freestream const &fs);
+ // Adjoint
+ static Eigen::VectorXd buildGradientAoa(tbox::Element const &e, std::vector<double> const &phi, Freestream const &fs);
+};
+
+} // namespace flow
+#endif //WFREESTREAMRESIDUAL_H
diff --git a/flow/src/wKuttaElement.cpp b/flow/src/wKuttaElement.cpp
index 16a9f5fe..798fd778 100644
--- a/flow/src/wKuttaElement.cpp
+++ b/flow/src/wKuttaElement.cpp
@@ -51,170 +51,6 @@ KuttaElement::KuttaElement(size_t _no, Element *&_surE, Element *&_volE,
nRow = teN.size();
}
-/**
- * @brief Compute stiffness matrices associated to Kutta
- *
- * K_ij = int V*dN_j * N_i dS
- *
- * 3
- * o
- * / \
- * / \
- * / \
- * 1,1' o-x---x-o 2,2'
- * A B
- *
- * I = sum_j w'_j psi_j N'_j V_j invJ_j dN_j phi_j dtm(J')_j
- * = psi sum_j (w'_j N'_j V_j invJ_j dN_j dtm(J')_j) phi
- * = (psi)|1X2 (sum_j K_j)|2X3 (phi)|3X1
- *
- * K = w'[N'1][Vx Vy]inv([dx/dxi dy/dxi])[dxiN1 dxiN2 dxiN3] dtm([dx/dxi' dy/dxi']^2)
- * [N'2] ([dx/det dy/det])[detN1 detN2 detN3]
- *
- * phi1 phi2 phi3
- * K = 1' [ . . . ]
- * 2' [ . . . ]
- *
- * NB: K should be eval at GP A and B, for upper and lower element
- * NB: V = grad(phi)
- * @authors Adrien Crovato
- */
-void KuttaElement::buildK(Eigen::MatrixXd &K, std::vector<double> const &u)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- Eigen::MatrixXd const &volDff = volE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &surMem = surE->getVMem();
- Eigen::MatrixXd const &volJ = volE->getVMem().getJinv(0);
- // Reset matrix
- K = Eigen::MatrixXd::Zero(nRow, nCol);
-
- // Volume contribution
- // V & dN
- // A = V*[inv(J)*dN]
- Eigen::RowVectorXd A(nCol);
- A = volE->computeGrad(u, 0).transpose() * volJ * volDff;
-
- // Build
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // N
- Eigen::VectorXd const &ff = surCache.getSf(k);
- Eigen::VectorXd N(nRow);
- for (size_t i = 0; i < nRow; ++i)
- N(i) = ff(teN[i].first);
-
- // K = K + N*B*w*dtm
- K += N * A * sign * surGauss.getW(k) * surMem.getDetJ(k);
- }
-}
-
-/**
- * @brief Compute stiffness matrices associated to Kutta for Newton-Raphson
- *
- * K_ij = int 2*V*dN_j * N_i dS
- * @authors Adrien Crovato
- */
-void KuttaElement::buildNK(Eigen::MatrixXd &K, std::vector<double> const &u)
-{
- // Build
- buildK(K, u);
- K *= 2.;
-}
-
-/**
- * @brief Compute flux vectors associated to Kutta for Newton-Raphson
- *
- * s_i = int V*V * N_i dS
- * @authors Adrien Crovato
- */
-void KuttaElement::buildNs(Eigen::VectorXd &b, std::vector<double> const &u)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- // Get Jacobian
- Mem &surMem = surE->getVMem();
- // Reset matrices
- b = Eigen::VectorXd::Zero(nRow);
-
- // Build
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // N
- Eigen::VectorXd const &ff = surCache.getSf(k);
- Eigen::VectorXd N(nRow);
- for (size_t i = 0; i < nRow; ++i)
- N(i) = ff(teN[i].first);
- // V2
- double grad2 = volE->computeGrad2(u, 0);
-
- // Compute flux vector
- b += sign * N * grad2 * surGauss.getW(k) * surMem.getDetJ(k);
- }
-}
-
-/**
- * @brief Compute stiffness matrices associated to Kutta for Newton-Raphson
- *
- * K_ij = int d(2*V*dN_j * N_i) dS
- * = sum 2*V*dgrad_phi * N_i * detJ
- * + sum grad_phi^2 * N_i ddetJ
- */
-void KuttaElement::buildGradientResidualMesh(Eigen::MatrixXd &K, Eigen::MatrixXd &Ks, std::vector<double> const &u)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- // Get Jacobian
- Mem &surMem = surE->getVMem();
- Mem &volMem = volE->getVMem();
- Eigen::MatrixXd const &iJ = volMem.getJinv(0);
- std::vector<Eigen::MatrixXd> const &dJ = volMem.getGradJ(0);
- // Reset matrices
- K = Eigen::MatrixXd::Zero(surE->nodes.size(), nDim * volE->nodes.size());
- Ks = Eigen::MatrixXd::Zero(surE->nodes.size(), nDim * surE->nodes.size());
-
- // Gradient of potential
- Eigen::VectorXd dPhi = volE->computeGrad(u, 0);
-
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // weight and Jacobian determinant
- double w = surGauss.getW(k);
- double detJ = surMem.getDetJ(k);
- std::vector<double> dDetJ;
- surE->buildGradientJ(k, dDetJ);
- // shape functions
- Eigen::VectorXd const &sf = surCache.getSf(k);
- Eigen::VectorXd n(nRow);
- for (size_t i = 0; i < nRow; ++i)
- n(i) = sf(teN[i].first);
-
- for (size_t i = 0; i < volE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- // derivative of velocity jump
- K.col(idx) += sign * w * n * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // psi * dgrad_phi^2
- }
- }
- // Surface
- for (size_t i = 0; i < surE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- // derivative of jacobian determinant
- Ks.col(idx) += sign * w * n * dPhi.dot(dPhi) * dDetJ[idx]; // psi * grad_phi^2 * ddetJ
- }
- }
- }
-}
-
void KuttaElement::write(std::ostream &out) const
{
out << "KuttaElement #" << no
diff --git a/flow/src/wKuttaElement.h b/flow/src/wKuttaElement.h
index 9a15f94a..4103cc39 100644
--- a/flow/src/wKuttaElement.h
+++ b/flow/src/wKuttaElement.h
@@ -48,23 +48,18 @@ public:
tbox::Element *surE; ///< Surface element
tbox::Element *volE; ///< Connected volume element
- size_t nRow; ///< number of rows in stiffness matrix
- size_t nCol; ///< number of columns in stiffness matrix
+ size_t nRow; ///< number of rows in stiffness matrix
+ size_t nCol; ///< number of columns in stiffness matrix
+ size_t nDim; ///< dimension of volume element
+
std::vector<std::pair<size_t, tbox::Node *>> teN; ///< Map of local node indices and trailing edge nodes
+ int sign; ///< +1 if upper element, -1 if lower element
KuttaElement(size_t _no, tbox::Element *&_surE, tbox::Element *&_volE, std::vector<std::pair<size_t, tbox::Node *>> &_teN, int _sign);
#ifndef SWIG
- void buildK(Eigen::MatrixXd &K, std::vector<double> const &u);
- void buildNK(Eigen::MatrixXd &K, std::vector<double> const &u);
- void buildNs(Eigen::VectorXd &b, std::vector<double> const &u);
- void buildGradientResidualMesh(Eigen::MatrixXd &K, Eigen::MatrixXd &Ks, std::vector<double> const &u);
virtual void write(std::ostream &out) const override;
#endif
-
-private:
- size_t nDim; ///< dimension of volume element
- int sign; ///< +1 if upper element, -1 if lower element
};
} // namespace flow
diff --git a/flow/src/wKuttaResidual.cpp b/flow/src/wKuttaResidual.cpp
new file mode 100644
index 00000000..59c2010a
--- /dev/null
+++ b/flow/src/wKuttaResidual.cpp
@@ -0,0 +1,167 @@
+/*
+ * 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 "wKuttaResidual.h"
+#include "wKuttaElement.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Build the residual matrix for a fixed-point iteration, on one Kutta element
+ *
+ * A = \int psi * v_u*grad_phi dS
+ */
+Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<double> const &phi)
+{
+ // Get shape functions, Jacobian and Gauss points
+ Cache &sCache = ke.surE->getVCache();
+ Gauss &sGauss = sCache.getVGauss();
+ Mem &sMem = ke.surE->getVMem();
+ // Get size
+ size_t nRow = ke.nRow;
+
+ // Velocity
+ Eigen::RowVectorXd Aup = ke.volE->computeGrad(phi, 0).transpose() * ke.volE->getVMem().getJinv(0) * ke.volE->getVCache().getDsf(0);
+
+ // Build
+ Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
+ for (size_t k = 0; k < sGauss.getN(); ++k)
+ {
+ // shape functions
+ Eigen::VectorXd const &sf = sCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = sf(ke.teN[i].first);
+ // wake contribution
+ A += ke.sign * sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
+ }
+ return A;
+}
+
+/**
+ * @brief Build the residual vector, on one Kutta element
+ *
+ * b = \int psi * (grad_phi)^2 dS
+ */
+Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> const &phi)
+{
+ // Get shape functions, Jacobian and Gauss points
+ Cache &sCache = ke.surE->getVCache();
+ Gauss &sGauss = sCache.getVGauss();
+ Mem &sMem = ke.surE->getVMem();
+ // Get size
+ size_t nRow = ke.nRow;
+
+ // Velocity
+ double dPhi2 = ke.volE->computeGrad2(phi, 0);
+
+ // Build
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(nRow);
+ for (size_t k = 0; k < sGauss.getN(); ++k)
+ {
+ // shape functions
+ Eigen::VectorXd const &sf = sCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = sf(ke.teN[i].first);
+ // wake contribution
+ b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * sMem.getDetJ(k);
+ }
+ return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one Kutta element
+ *
+ * A = \int psi * 2*(grad_phi_u)*dgrad_phi_u dS
+ */
+Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::vector<double> const &phi)
+{
+ // Build
+ Eigen::MatrixXd A = KuttaResidual::buildFixed(ke, phi);
+ return 2 * A;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the nodes, on one Kutta element
+ *
+ * A = d( \int psi * (grad_phi)^2 dS )
+ * = \int psi * d(grad_phi)^2 dS
+ * + \int psi * (grad_phi)^2 ddS
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(KuttaElement const &ke, std::vector<double> const &phi)
+{
+ // Get shape functions and Gauss points
+ Cache &surCache = ke.surE->getVCache();
+ Gauss &surGauss = surCache.getVGauss();
+ // Get Jacobian
+ Mem &surMem = ke.surE->getVMem();
+ Mem &volMem = ke.volE->getVMem();
+ Eigen::MatrixXd const &iJ = volMem.getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJ = volMem.getGradJ(0);
+ // Get size
+ size_t nRow = ke.nRow;
+ size_t nDim = ke.nDim;
+ size_t nCol = ke.nCol;
+ size_t nSur = ke.surE->nodes.size();
+
+ // velocity
+ Eigen::VectorXd dPhi = ke.volE->computeGrad(phi, 0);
+
+ // Build
+ Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(nRow, nDim * nCol);
+ Eigen::MatrixXd As = Eigen::MatrixXd::Zero(nRow, nDim * nSur);
+ for (size_t k = 0; k < surGauss.getN(); ++k)
+ {
+ // weight and Jacobian determinant
+ double w = surGauss.getW(k);
+ double detJ = surMem.getDetJ(k);
+ std::vector<double> dDetJ = surMem.getGradDetJ(k);
+ // shape functions
+ Eigen::VectorXd const &sf = surCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = sf(ke.teN[i].first);
+
+ // Volume
+ for (size_t i = 0; i < nCol; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ // gradient of velocity jump
+ Av.col(idx) += ke.sign * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // psi * dgrad_phi^2
+ }
+ }
+ // Surface
+ for (size_t i = 0; i < nSur; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ // gradient of jacobian determinant
+ As.col(idx) += ke.sign * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // psi * grad_phi^2 * ddetJ
+ }
+ }
+ }
+ return std::make_tuple(Av, As);
+}
diff --git a/flow/src/wKuttaResidual.h b/flow/src/wKuttaResidual.h
new file mode 100644
index 00000000..562b3dcc
--- /dev/null
+++ b/flow/src/wKuttaResidual.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.
+ */
+
+#ifndef WKUTTARESIDUAL_H
+#define WKUTTARESIDUAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of Kutta residuals
+ */
+class FLOW_API KuttaResidual
+{
+public:
+ // Picard
+ static Eigen::MatrixXd buildFixed(KuttaElement const &ke, std::vector<double> const &phi);
+ // Newton
+ static Eigen::VectorXd build(KuttaElement const &ke, std::vector<double> const &phi);
+ static Eigen::MatrixXd buildGradientFlow(KuttaElement const &ke, std::vector<double> const &phi);
+ // Adjoint
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(KuttaElement const &ke, std::vector<double> const &phi);
+};
+
+} // namespace flow
+#endif //WKUTTARESIDUAL_H
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
new file mode 100644
index 00000000..a251f54f
--- /dev/null
+++ b/flow/src/wLoadFunctional.cpp
@@ -0,0 +1,136 @@
+/*
+ * 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 "wLoadFunctional.h"
+#include "wFct0.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Compute the loads due to pressure normalized by the freestream dynamic pressure, on one boundary element
+ *
+ * b = \int psi * cp * n dS
+ */
+Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Pressure force coefficient in associated volume
+ Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.normal();
+
+ // Integrate cp on the surface
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(3 * e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ Eigen::VectorXd sf = cache.getSf(k);
+ double wdj = gauss.getW(k) * mem.getDetJ(k);
+ for (size_t i = 0; i < e.nodes.size(); ++i)
+ b.block<3, 1>(i * 3, 0) += sf(i) * fcp * wdj; // stack nodal force vectors [fx1, fy1, fz1, fx2, fy2, fz2, ...]
+ }
+ return b;
+}
+
+/**
+ * @brief Build the gradient of the loads with respect to the flow variable, on one boundary element
+ *
+ * A = \int psi * dcp * n dS
+ */
+Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Gradient of pressure force coefficient in associated volume
+ Eigen::MatrixXd dfcp = e.normal() * cp.evalD(eV, phi, 0) * eV.computeGrad(phi, 0).transpose() * eV.getVMem().getJinv(0) * eV.getVCache().getDsf(0);
+
+ // Integrate cp on the surface
+ Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ Eigen::VectorXd sf = cache.getSf(k);
+ double wdj = gauss.getW(k) * mem.getDetJ(k);
+ for (size_t i = 0; i < e.nodes.size(); ++i)
+ A.block(i * 3, 0, 3, eV.nodes.size()) += sf(i) * dfcp * wdj;
+ }
+ return A;
+}
+
+/**
+ * @brief Build the gradient of the loads with respect to the nodes, on one boundary element
+ *
+ * A = d( \int psi * cp * n dS )
+ * = \int psi * dcp * n * dS
+ * + \int psi * cp * dn * dS
+ * + \int psi * cp * n * ddS
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp, int nDim)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Variables
+ double cP = cp.eval(eV, phi, 0); // pressure coefficient
+ Eigen::Vector3d nrm = e.normal(); // unit normal
+ // Gradients
+ double dCp = cp.evalD(eV, phi, 0); // pressure coefficient
+ Eigen::VectorXd dPhi = eV.computeGrad(phi, 0); // potential
+ Eigen::MatrixXd const &iJ = eV.getVMem().getJinv(0); // inverse Jacobian
+ std::vector<Eigen::MatrixXd> const &dJ = eV.getVMem().getGradJ(0); // Jacobian
+ std::vector<Eigen::Vector3d> dNrm = e.gradientNormal(); // unit normal
+
+ // Build element contributions
+ Eigen::MatrixXd As = Eigen::MatrixXd::Zero(3 * e.nodes.size(), e.nodes.size() * nDim);
+ Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size() * nDim);
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ Eigen::VectorXd sf = cache.getSf(k);
+ double w = gauss.getW(k);
+ double detJ = mem.getDetJ(k);
+ std::vector<double> dDetJ = mem.getGradDetJ(k);
+ for (size_t i = 0; i < e.nodes.size(); ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ // surface
+ for (size_t j = 0; j < e.nodes.size(); ++j)
+ {
+ size_t idx = j * nDim + m;
+ As.block<3, 1>(i * 3, idx) += sf(i) * cP * w * (dNrm[idx] * detJ + nrm * dDetJ[idx]); // Cp * dn * S + Cp * n * dS
+ }
+ // volume
+ for (size_t j = 0; j < eV.nodes.size(); ++j)
+ {
+ size_t idx = j * nDim + m;
+ Av.block<3, 1>(i * 3, idx) += sf(i) * dCp * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * nrm * w * detJ; // dCp * n * S
+ }
+ }
+ }
+ }
+ return std::make_tuple(As, Av);
+}
diff --git a/flow/src/wLoadFunctional.h b/flow/src/wLoadFunctional.h
new file mode 100644
index 00000000..4a1a25a8
--- /dev/null
+++ b/flow/src/wLoadFunctional.h
@@ -0,0 +1,42 @@
+/*
+ * 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 WLOADFUNCTIONAL_H
+#define WLOADFUNCTIONAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of load functional
+ */
+class FLOW_API LoadFunctional
+{
+public:
+ // Solver
+ static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp);
+ // Adjoint
+ static Eigen::MatrixXd buildGradientFlow(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, tbox::Fct0 const &cp, int nDim);
+};
+
+} // namespace flow
+#endif //WLOADFUNCTIONAL_H
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index 4b24dc08..4ac08c51 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -23,13 +23,17 @@
#include "wKutta.h"
#include "wBlowing.h"
#include "wFpeLSFunction.h"
+#include "wPotentialResidual.h"
+#include "wFreestreamResidual.h"
+#include "wWakeResidual.h"
+#include "wKuttaResidual.h"
+#include "wBlowingResidual.h"
#include "wMshData.h"
#include "wTag.h"
#include "wNode.h"
#include "wElement.h"
#include "wMem.h"
-#include "wCache.h"
#include "wLinesearch.h"
#include "wLinearSolver.h"
@@ -66,7 +70,7 @@ Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
/**
* @brief Run the Newton solver
- *
+ *
* Solve the steady transonic Full Potential Equation
* @authors Adrien Crovato
*/
@@ -268,60 +272,35 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
Element *e = p.first;
//Upwind element
Element *eU = p.second[0];
+ // Build elementary matrices
+ Eigen::MatrixXd Ae1, Ae2;
+ std::tie(Ae1, Ae2) = PotentialResidual::buildGradientFlow(*e, *eU, phi, *fluid, muCv, mCOv);
- // Subsonic contribution (drho*grad_phi*grad_psi + rho*grad_phi*grad_psi)
- Eigen::MatrixXd Ae1 = e->buildDK(phi, fluid->rho);
- Eigen::MatrixXd Ae2 = e->buildK(phi, fluid->rho);
-
- // Supersonic contribution
+ // Assembly (subsonic)
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae1(ii, jj)));
+ }
+ }
+ // Assembly (supersonic)
double mach = fluid->mach.eval(*e, phi, 0);
if (mach > mCOv)
{
- // switching function and derivative
- double mu = muCv * (1 - (mCOv * mCOv) / (mach * mach));
- double dmu = 2 * muCv * mCOv * mCOv / (mach * mach * mach);
-
- // scale 1st and 2nd terms
- Ae1 *= 1 - mu;
- Ae2 *= 1 - mu;
-
- // 3rd term (mu*drhoU*grad_phi*grad_psi)
- Eigen::MatrixXd Ae3 = mu * fluid->rho.evalD(*eU, phi, 0) * e->buildDs(phi, Fct0C(1.)) * eU->computeGrad(phi, 0).transpose() * eU->getVMem().getJinv(0) * eU->getVCache().getDsf(0);
-
- // 4th term (mu*rhoU*grad_phi*grad_psi)
- Eigen::MatrixXd Ae4 = mu * fluid->rho.eval(*eU, phi, 0) * e->buildK(phi, Fct0C(1.));
-
- // 5th term dmu*(rhoU-rho)*grad_phi*grad_psi
- Eigen::MatrixXd Ae5 = dmu * (fluid->rho.eval(*eU, phi, 0) - fluid->rho.eval(*e, phi, 0)) * e->buildDK(phi, fluid->mach);
-
- // Assembly (supersonic)
- tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
- {
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae4(ii, jj) + Ae5(ii, jj)));
- }
for (size_t jj = 0; jj < eU->nodes.size(); ++jj)
{
Node *nodj = eU->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae3(ii, jj)));
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae2(ii, jj)));
}
}
}
- // Assembly (subsonic)
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
- {
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ae1(ii, jj) + Ae2(ii, jj)));
- }
- }
});
tms["0-Jbase"].stop();
@@ -335,7 +314,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build Ku, Kl
Eigen::MatrixXd Aue, Ale;
- we->buildNK(Aue, Ale, phi, vi);
+ std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi, vi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -361,8 +340,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
{
tbb::parallel_do(kutta->kEle.begin(), kutta->kEle.end(), [&](KuttaElement *ke) {
// Build K
- Eigen::MatrixXd Ae;
- ke->buildNK(Ae, phi);
+ Eigen::MatrixXd Ae = KuttaResidual::buildGradientFlow(*ke, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < ke->nRow; ++ii)
@@ -425,28 +403,8 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
Element *e = p.first;
// Upwind element
Element *eU = p.second[0];
-
- // Subsonic contribution
- Eigen::VectorXd be = e->buildDs(phi, fluid->rho);
-
- // Supersonic contribution
- double mach = fluid->mach.eval(*e, phi, 0);
- if (mach > mCOv)
- {
- double mu = muCv * (1 - (mCOv * mCOv) / (mach * mach));
- double rhoU = fluid->rho.eval(*eU, phi, 0);
- Eigen::VectorXd beU = e->buildDs(phi, Fct0C(1.));
- // scale 1st and 2nd terms
- be *= 1 - mu;
- beU *= mu * rhoU;
- // Assembly (supersonic)
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- R(rows[nodi->row]) -= beU(ii);
- }
- }
+ // Build elementary vector
+ Eigen::VectorXd be = PotentialResidual::build(*e, *eU, phi, *fluid, muCv, mCOv);
// Assembly (subsonic)
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
@@ -460,7 +418,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
{
tbb::parallel_do(nBC->tag->elems.begin(), nBC->tag->elems.end(), [&](Element *e) {
// Build elementary flux vector
- Eigen::VectorXd be = e->builds(phi, nBC->f);
+ Eigen::VectorXd be = FreestreamResidual::build(*e, phi, *nBC);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -476,7 +434,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
{
tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) {
// Build elementary flux vector
- Eigen::VectorXd be = p.first->builds(phi, *(p.second));
+ Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -498,7 +456,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build bu, bl
Eigen::VectorXd bue, ble;
- we->buildNs(bue, ble, phi, vi);
+ std::tie(bue, ble) = WakeResidual::build(*we, phi, vi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -516,8 +474,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
{
tbb::parallel_do(kutta->kEle.begin(), kutta->kEle.end(), [&](KuttaElement *ke) {
// Build b
- Eigen::VectorXd be;
- ke->buildNs(be, phi);
+ Eigen::VectorXd be = KuttaResidual::build(*ke, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < ke->nRow; ++ii)
diff --git a/flow/src/wPicard.cpp b/flow/src/wPicard.cpp
index 7849d56c..8b66292b 100644
--- a/flow/src/wPicard.cpp
+++ b/flow/src/wPicard.cpp
@@ -22,6 +22,11 @@
#include "wWake.h"
#include "wKutta.h"
#include "wBlowing.h"
+#include "wPotentialResidual.h"
+#include "wFreestreamResidual.h"
+#include "wWakeResidual.h"
+#include "wKuttaResidual.h"
+#include "wBlowingResidual.h"
#include "wMshData.h"
#include "wTag.h"
@@ -58,7 +63,7 @@ Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
/**
* @brief Run the Picard sovler
- *
+ *
* Solve the steady subsonic Full Potential Equation
* @authors Adrien Crovato
*/
@@ -199,7 +204,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
Element *e = p.first;
// Elementary stiffness
- Eigen::MatrixXd Ae = e->buildK(phi, fluid->rho);
+ Eigen::MatrixXd Ae = PotentialResidual::buildFixed(*e, phi, *fluid);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -225,7 +230,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build Ku, Kl
Eigen::MatrixXd Aue, Ale;
- we->buildK(Aue, Ale, phi, vi);
+ std::tie(Aue, Ale) = WakeResidual::buildFixed(*we, phi, vi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -251,8 +256,8 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
{
tbb::parallel_do(kutta->kEle.begin(), kutta->kEle.end(), [&](KuttaElement *ke) {
// Build K
- Eigen::MatrixXd Ae;
- ke->buildK(Ae, phi);
+ Eigen::MatrixXd Ae = KuttaResidual::buildFixed(*ke, phi);
+
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < ke->nRow; ++ii)
@@ -274,7 +279,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
{
tbb::parallel_do(nBC->tag->elems.begin(), nBC->tag->elems.end(), [&](Element *e) {
// Build elementary flux vector
- Eigen::VectorXd be = e->builds(phi, nBC->f);
+ Eigen::VectorXd be = FreestreamResidual::build(*e, phi, *nBC);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -290,7 +295,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
{
tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) {
// Build elementary flux vector
- Eigen::VectorXd be = p.first->builds(phi, *(p.second));
+ Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp
new file mode 100644
index 00000000..ff205478
--- /dev/null
+++ b/flow/src/wPotentialResidual.cpp
@@ -0,0 +1,250 @@
+/*
+ * 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 "wPotentialResidual.h"
+#include "wMedium.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Build the residual matrix for a fixed-point iteration, on one element
+ *
+ * A = \int rho * grad_phi * grad_psi dV
+ */
+Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<double> const &phi, Medium const &fluid)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Build
+ Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Shape functions gradient
+ Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ // Elementary stiffness matrix
+ A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return A;
+}
+
+/**
+ * @brief Build the residual vector, on one element
+ *
+ * b = \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV
+ * @todo remove muC, mCO + mu as function?
+ */
+Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Subsonic contribution
+ Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ b += fluid.rho.eval(e, phi, k) * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGrad(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+
+ // Supersonic contribution
+ double mach = fluid.mach.eval(e, phi, 0);
+ if (mach > mCO)
+ {
+ double mu = muC * (1 - (mCO * mCO) / (mach * mach)); // switching function
+ double rhoU = fluid.rho.eval(eU, phi, 0); // density on upwind element
+ // scale 1st term
+ b *= 1 - mu;
+ // contribution of upwind element
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ b += mu * rhoU * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGrad(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+ }
+ return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one element
+ *
+ * A = \int d( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV
+ * = \int (1-mu)*drho * grad_phi * grad_psi dV
+ * + \int mu*drhoU * grad_phi * grad_psi dV
+ * + \int ( (1-mu)*rho + mu*rhoU ) * dgrad_phi * grad_psi dV
+ * + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV
+ * @todo remove muC, mCO + mu as function?
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Subsonic contribution
+ Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Gauss point and determinant
+ double wdj = gauss.getW(k) * mem.getDetJ(k);
+ // Shape functions and solution gradient
+ Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+
+ // rho * grad_phi*grad_psi + drho * grad_phi*grad_psi
+ A1 += fluid.rho.eval(e, phi, k) * dSf.transpose() * dSf * wdj;
+ A1 += fluid.rho.evalD(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
+ }
+
+ // Supersonic contribution
+ Eigen::MatrixXd A2;
+ double mach = fluid.mach.eval(e, phi, 0);
+ if (mach > mCO)
+ {
+ // switching function and gradient
+ double mu = muC * (1 - (mCO * mCO) / (mach * mach));
+ double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
+ // density and gradients on upwind elements
+ Eigen::VectorXd dPhiU = eU.computeGrad(phi, 0);
+ double rhoU = fluid.rho.eval(eU, phi, 0);
+ double dRhoU = fluid.rho.evalD(eU, phi, 0);
+ Eigen::MatrixXd const &dSfU = eU.getVMem().getJinv(0) * eU.getVCache().getDsf(0);
+
+ // upwinding
+ A1 *= 1 - mu;
+ A2 = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Gauss point and determinant
+ double wdj = gauss.getW(k) * mem.getDetJ(k);
+ // Shape functions and solution gradient
+ Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+ Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+
+ // mu*drhoU*grad_phi*grad_psi
+ A2 += mu * dRhoU * dSf.transpose() * dPhi * dPhiU.transpose() * dSfU * wdj;
+ // mu*rhoU*grad_phi*grad_psi + dmu*(rhoU-rho)*grad_phi*grad_psi
+ A1 += mu * rhoU * dSf.transpose() * dSf * wdj;
+ A1 += dmu * (rhoU - fluid.rho.eval(e, phi, k)) * fluid.mach.evalD(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
+ }
+ }
+ return std::make_tuple(A1, A2);
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the nodes, on one element
+ *
+ * A = d( \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV )
+ * = \int (1-mu)*drho * grad_phi * grad_psi dV
+ * + \int mu*drhoU * grad_phi * grad_psi dV
+ * + \int ( (1-mu)*rho + mu*rhoU ) * dgrad_phi * grad_psi dV
+ * + \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * dgrad_psi dV
+ * + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV
+ * + \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi ddV
+ * @todo remove muC, mCO + mu as function?
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, int nDim, double muC, double mCO)
+{
+ // Get pre-computed values
+ Cache &cache = e.getVCache();
+ Gauss &gauss = cache.getVGauss();
+ Mem &mem = e.getVMem();
+
+ // Subsonic contribution
+ Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Density and gradients
+ Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ double rho = fluid.rho.eval(e, phi, k);
+ double dRho = fluid.rho.evalD(e, phi, k);
+ Eigen::MatrixXd const &dSf = cache.getDsf(k);
+ Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ // Jacobian gradients
+ std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
+ std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+ // Gauss points and determinant
+ double w = gauss.getW(k);
+ double detJ = mem.getDetJ(k);
+ for (size_t i = 0; i < e.nodes.size(); ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ A1.col(idx) += w * dRho * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // drho * grad_phi*grad_psi * detj
+ A1.col(idx) += w * rho * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // rho * dgrad_phi * grad_psi * detj
+ A1.col(idx) += w * rho * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // rho * grad_phi * dgrad_psi * detj
+ A1.col(idx) += w * rho * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // rho * grad_phi * grad_psi * ddetj
+ }
+ }
+ }
+ // Supersonic contribution
+ Eigen::MatrixXd A2;
+ double mach = fluid.mach.eval(e, phi, 0);
+ if (mach > mCO)
+ {
+ // switching function and gradient
+ double mu = muC * (1 - (mCO * mCO) / (mach * mach));
+ double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
+ // upwind density and gradient
+ Eigen::VectorXd dPhiU = eU.computeGrad(phi, 0);
+ double rhoU = fluid.rho.eval(eU, phi, 0);
+ double dRhoU = fluid.rho.evalD(eU, phi, 0);
+ Eigen::MatrixXd const &iJU = eU.getVMem().getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJU = eU.getVMem().getGradJ(0);
+
+ // upwinding
+ A1 *= 1 - mu;
+ A2 = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * eU.nodes.size());
+ for (size_t k = 0; k < gauss.getN(); ++k)
+ {
+ // Density, Mach and Gradients
+ Eigen::VectorXd dPhi = e.computeGrad(phi, k);
+ double rho = fluid.rho.eval(e, phi, k);
+ double dM = fluid.mach.evalD(e, phi, k);
+ Eigen::MatrixXd const &dSf = cache.getDsf(k);
+ Eigen::MatrixXd const &iJ = mem.getJinv(k);
+ // Jacobian gradients
+ std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
+ std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+ // Gauss points and determinant
+ double w = gauss.getW(k);
+ double detJ = mem.getDetJ(k);
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ for (size_t i = 0; i < e.nodes.size(); ++i)
+ {
+ size_t idx = i * nDim + m;
+ A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // mu * rhoU * dgrad_phi * grad_psi * detj
+ A1.col(idx) += w * mu * rhoU * (-iJ * dJ[idx] * iJ * dSf).transpose() * dPhi * detJ; // mu * rhoU * grad_phi * dgrad_psi * detj
+ A1.col(idx) += w * mu * rhoU * (iJ * dSf).transpose() * dPhi * dDetJ[idx]; // mu * rhoU * grad_phi * grad_psi * ddetj
+ A1.col(idx) += w * (-rho + rhoU) * dmu * dM * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * (iJ * dSf).transpose() * dPhi * detJ; // (rhoU-rho) * dmu * grad_phi * grad_psi * detj
+ }
+ for (size_t i = 0; i < eU.nodes.size(); ++i)
+ {
+ size_t idx = i * nDim + m;
+ A2.col(idx) += w * mu * dRhoU * dPhiU.transpose() * (-iJU * dJU[idx] * dPhiU) * (iJ * dSf).transpose() * dPhi * detJ; // mu * drhoU * grad_phi * grad_psi * detj
+ }
+ }
+ }
+ }
+ return std::make_tuple(A1, A2);
+}
diff --git a/flow/src/wPotentialResidual.h b/flow/src/wPotentialResidual.h
new file mode 100644
index 00000000..62dde67b
--- /dev/null
+++ b/flow/src/wPotentialResidual.h
@@ -0,0 +1,45 @@
+/*
+ * 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 WPOTENTIALRESIDUAL_H
+#define WPOTENTIALRESIDUAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of nonlinear potential equation residuals
+ * @todo split stabilization terms from subsonic terms?
+ */
+class FLOW_API PotentialResidual
+{
+public:
+ // Picard
+ static Eigen::MatrixXd buildFixed(tbox::Element const &e, std::vector<double> const &phi, Medium const &fluid);
+ // Newton
+ static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO);
+ // Adjoint
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, int nDim, double muC, double mCO);
+};
+
+} // namespace flow
+#endif //WPOTENTIALRESIDUAL_H
diff --git a/flow/src/wProblem.cpp b/flow/src/wProblem.cpp
index 2babefbc..966064ca 100644
--- a/flow/src/wProblem.cpp
+++ b/flow/src/wProblem.cpp
@@ -38,6 +38,55 @@ Problem::Problem(std::shared_ptr<MshData> _msh, int dim, double aoa, double aos,
x_ref(2) = zref;
}
+/**
+ * @brief Compute drag, sideforce and lift normalized directions
+ * @todo use F1ElVi
+ */
+std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> Problem::computeDirections()
+{
+ Eigen::Vector3d dirD, dirS, dirL;
+ if (nDim == 2)
+ {
+ dirD = Eigen::Vector3d(cos(alpha), sin(alpha), 0) / S_ref;
+ dirS = Eigen::Vector3d(0, 0, 0);
+ dirL = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / S_ref;
+ }
+ else
+ {
+ dirD = Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)) / S_ref;
+ dirS = Eigen::Vector3d(-cos(alpha) * sin(beta), cos(beta), -sin(alpha) * sin(beta)) / S_ref;
+ dirL = Eigen::Vector3d(-sin(alpha), 0, cos(alpha)) / S_ref;
+ }
+ return std::make_tuple(dirD, dirS, dirL);
+}
+
+/**
+ * @brief Compute gradient of drag, sideforce and lift normalized directions
+ * @todo use F1ElVi
+ * @todo only wrt alpha
+ */
+std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> Problem::computeGradientDirections()
+{
+ Eigen::Vector3d dirD, dirS, dirL;
+ if (nDim == 2)
+ {
+ dirD = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / S_ref;
+ dirS = Eigen::Vector3d(0, 0, 0);
+ dirL = Eigen::Vector3d(-cos(alpha), -sin(alpha), 0) / S_ref;
+ }
+ else
+ {
+ dirD = Eigen::Vector3d(-sin(alpha) * cos(beta), sin(beta), cos(alpha) * cos(beta)) / S_ref;
+ dirS = Eigen::Vector3d(sin(alpha) * sin(beta), cos(beta), -cos(alpha) * sin(beta)) / S_ref;
+ dirL = Eigen::Vector3d(-cos(alpha), 0, -sin(alpha)) / S_ref;
+ }
+ return std::make_tuple(dirD, dirS, dirL);
+}
+
+/**
+ * @brief Compute freestream velocity
+ * @todo use F1ElVi
+ */
Eigen::VectorXd Problem::computeVi()
{
if (nDim == 2)
@@ -46,6 +95,19 @@ Eigen::VectorXd Problem::computeVi()
return Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta));
}
+/**
+ * @brief Compute freestream velocity
+ * @todo use F1ElVi
+ * @todo only wrt alpha
+ */
+Eigen::VectorXd Problem::computeGradientVi()
+{
+ if (nDim == 2)
+ return Eigen::Vector2d(-sin(alpha), cos(alpha));
+ else
+ return Eigen::Vector3d(-sin(alpha) * cos(beta), 0, cos(alpha) * cos(beta));
+}
+
/**
* @brief Set the fluid medium
* @authors Adrien Crovato
diff --git a/flow/src/wProblem.h b/flow/src/wProblem.h
index 008db76e..3f4b9b63 100644
--- a/flow/src/wProblem.h
+++ b/flow/src/wProblem.h
@@ -62,7 +62,10 @@ public:
double zref);
virtual ~Problem() { std::cout << "~Problem()\n"; }
+ std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> computeDirections();
+ std::tuple<Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d> computeGradientDirections();
Eigen::VectorXd computeVi();
+ Eigen::VectorXd computeGradientVi();
void set(std::shared_ptr<Medium> m);
void add(std::shared_ptr<Boundary> b);
void set(std::shared_ptr<Initial> i);
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 6f7d494e..c90dda9f 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -20,11 +20,11 @@
#include "wAssign.h"
#include "wBoundary.h"
#include "wWake.h"
+#include "wLoadFunctional.h"
#include "wMshData.h"
#include "wNode.h"
#include "wElement.h"
-#include "wMem.h"
#include "wTag.h"
#include "wResults.h"
#include "wMshExport.h"
@@ -160,55 +160,50 @@ void Solver::save(MshExport *mshWriter, int n)
*/
void Solver::computeLoad()
{
+ // Reset coefficients
Cl = 0;
Cd = 0;
Cs = 0;
Cm = 0;
+
+ // Lift, drag and sideforce normalized directions
+ Eigen::Vector3d dirL, dirD, dirS;
+ std::tie(dirD, dirS, dirL) = pbl->computeDirections();
+
+ // Compute loads
for (auto sur : pbl->bnds)
{
- // Map each element to its value
- std::map<Element *, double> euMap;
- for (auto e : sur->groups[0]->tag->elems)
- euMap[e] = sur->integrate(e, phi, pbl->medium->cP);
+ // Reset
+ std::fill(sur->nLoads.begin(), sur->nLoads.end(), Eigen::Vector3d::Zero());
- // Compute normalized load (i.e. load / pressure)
- for (size_t i = 0; i < sur->nodes.size(); ++i)
+ // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element
+ for (auto e : sur->groups[0]->tag->elems)
{
- sur->cLoadX[i] = 0;
- sur->cLoadY[i] = 0;
- sur->cLoadZ[i] = 0;
- for (auto e : sur->neMap.at(sur->nodes[i]))
+ // Build nodal pressure load
+ Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, pbl->medium->cP);
+
+ // Assembly
+ for (size_t i = 0; i < e->nodes.size(); ++i)
{
- Eigen::Vector3d cLoad = euMap.at(e) * e->normal() / e->nodes.size();
- sur->cLoadX[i] += cLoad(0);
- sur->cLoadY[i] += cLoad(1);
- sur->cLoadZ[i] += cLoad(2);
+ size_t idx = sur->nMap.at(e->nodes[i]);
+ sur->nLoads[idx] += be.block<3, 1>(i * 3, 0);
}
}
-
- // 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)
+ // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
+ // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
sur->Cm = 0;
Eigen::Vector3d Cf(0, 0, 0);
- for (auto p : euMap)
+ for (size_t i = 0; i < sur->nodes.size(); ++i)
{
- Eigen::Vector3d l = p.first->getVMem().getCg() - pbl->x_ref; // lever arm
- Eigen::Vector3d cf = p.second * p.first->normal(); // force coefficient
+ Eigen::Vector3d l = sur->nodes[i]->pos - pbl->x_ref; // lever arm
+ Eigen::Vector3d cf = sur->nLoads[i]; // normalized force
Cf -= cf;
sur->Cm -= (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
}
// rotate to flow direction
- if (pbl->nDim == 2)
- {
- sur->Cd = (Cf(0) * cos(pbl->alpha) + Cf(1) * sin(pbl->alpha)) / pbl->S_ref;
- sur->Cl = (-Cf(0) * sin(pbl->alpha) + Cf(1) * cos(pbl->alpha)) / pbl->S_ref;
- }
- else
- {
- sur->Cd = (Cf(0) * cos(pbl->alpha) * cos(pbl->beta) + Cf(1) * sin(pbl->beta) + Cf(2) * sin(pbl->alpha) * cos(pbl->beta)) / pbl->S_ref;
- sur->Cs = (-Cf(0) * cos(pbl->alpha) * sin(pbl->beta) + Cf(1) * cos(pbl->beta) - Cf(2) * sin(pbl->alpha) * sin(pbl->beta)) / pbl->S_ref;
- sur->Cl = (-Cf(0) * sin(pbl->alpha) + Cf(2) * cos(pbl->alpha)) / pbl->S_ref;
- }
+ sur->Cd = Cf.dot(dirD);
+ sur->Cs = Cf.dot(dirS);
+ sur->Cl = Cf.dot(dirL);
// compute total
Cl += sur->Cl;
Cd += sur->Cd;
diff --git a/flow/src/wWakeElement.cpp b/flow/src/wWakeElement.cpp
index 6f946f9c..a75ccaa0 100644
--- a/flow/src/wWakeElement.cpp
+++ b/flow/src/wWakeElement.cpp
@@ -24,9 +24,9 @@ using namespace flow;
WakeElement::WakeElement(size_t _no, Element *&_surUpE, Element *&_surLwE,
Element *&_volUpE, Element *&_volLwE,
- std::vector<std::pair<size_t, Node *>> &_wkN) : wObject(), no(_no),
- surUpE(_surUpE), surLwE(_surLwE), volUpE(_volUpE), volLwE(_volLwE),
- wkN(_wkN)
+ std::vector<std::pair<size_t, Node *>> &_wkN)
+ : wObject(), no(_no), surUpE(_surUpE), surLwE(_surLwE),
+ volUpE(_volUpE), volLwE(_volLwE), wkN(_wkN)
{
// Sanity checks
if (surUpE->type() != surLwE->type())
@@ -70,7 +70,6 @@ WakeElement::WakeElement(size_t _no, Element *&_surUpE, Element *&_surLwE,
/**
* @brief Map the wake nodes to their local volume node indices
- * @authors Adrien Crovato
*/
void WakeElement::mapNodes()
{
@@ -95,292 +94,6 @@ void WakeElement::mapNodes()
}
}
-/**
- * @brief Compute stiffness matrices associated to Wake
- *
- * K_ij = int V*dN_j * V_inf*dN_i dS
- *
- * 3
- * o
- * / \
- * / \
- * / \
- * 1,1' o-x---x-o 2,2'
- * A B
- *
- * I = sum_j w'_j psi_j V_inf dN'_j V_j invJ_j dN_j phi_j dtm(J')_j
- * = psi sum_j (w'_j V_inf dN'_j V_j invJ_j dN_j dtm(J')_j) phi
- * = (psi)|1X2 (sum_j K_j)|2X3 (phi)|3X1
- *
- * K = w'[dxiN'1 detN'1]T(inv([dx/dxi dy/dxi]))[Vinf,x][Vx Vy]inv([dx/dxi dy/dxi])[dxiN1 dxiN2 dxiN3] dtm([dx/dxi' dy/dxi']^2)
- * [dxiN'1 detN'1] ( ([dx/det dy/det]))[Vinf,y] ([dx/det dy/det])[detN1 detN2 detN3]
- *
- * phi1 phi2 phi3
- * K = 1' [ . . . ]
- * 2' [ . . . ]
- *
- * NB: K should be eval at GP A and B, for upper and lower element
- * NB: V = grad(phi)
- * @authors Adrien Crovato
- */
-void WakeElement::buildK(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, std::vector<double> const &u, Eigen::VectorXd const &Vi)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surUpE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- Eigen::MatrixXd const &volUpDff = volUpE->getVCache().getDsf(0);
- Eigen::MatrixXd const &volLwDff = volLwE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &surMem = surUpE->getVMem();
- Eigen::MatrixXd const &volUpJ = volUpE->getVMem().getJinv(0);
- Eigen::MatrixXd const &volLwJ = volLwE->getVMem().getJinv(0);
- // Reset matrices
- Kup = Eigen::MatrixXd::Zero(nRow, nColUp);
- Klw = Eigen::MatrixXd::Zero(nRow, nColLw);
-
- // Volume contribution
- // dN' (stabilization)
- Eigen::MatrixXd dN(nRow, nDim);
- for (size_t i = 0; i < nRow; ++i)
- for (size_t j = 0; j < nDim; ++j)
- dN(i, j) = volUpDff(j, vsMap.at(wkN[i].first));
- Eigen::VectorXd JdNV(nRow);
- JdNV = dN * volUpJ.transpose() * Vi;
-
- // A = V*[inv(J)*dN]
- Eigen::RowVectorXd Aup(nColUp), Alw(nColLw);
- Aup = volUpE->computeGrad(u, 0).transpose() * volUpJ * volUpDff;
- Alw = volLwE->computeGrad(u, 0).transpose() * volLwJ * volLwDff;
-
- // Build
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // N
- Eigen::VectorXd const &ff = surCache.getSf(k);
- Eigen::VectorXd N(nRow);
- for (size_t i = 0; i < nRow; ++i)
- N(i) = 0.5 * sqrt(surMem.getVol()) * JdNV(i) + ff(wkN[i].first);
-
- // K = K + N*A*w*dtm
- Kup += N * Aup * surGauss.getW(k) * surMem.getDetJ(k);
- Klw += N * Alw * surGauss.getW(k) * surMem.getDetJ(k);
- }
-}
-
-/**
- * @brief Compute stiffness matrices associated to Kutta for Newton-Raphson
- *
- * K_ij = int 2*V*dN_j * V_inf*dN_i dS
- * @authors Adrien Crovato
- */
-void WakeElement::buildNK(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, std::vector<double> const &u, Eigen::VectorXd const &Vi)
-{
- // Build
- buildK(Kup, Klw, u, Vi);
- Kup *= 2.;
- Klw *= 2.;
-}
-
-/**
- * @brief Compute flux vectors associated to Kutta for Newton-Raphson
- *
- * s_i = int V*V * V_inf*dN_i dS
- * @authors Adrien Crovato
- */
-void WakeElement::buildNs(Eigen::VectorXd &bUp, Eigen::VectorXd &bLw, std::vector<double> const &u, Eigen::VectorXd const &Vi)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surUpE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- Eigen::MatrixXd const &volUpDff = volUpE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &surMem = surUpE->getVMem();
- Eigen::MatrixXd const &volUpJ = volUpE->getVMem().getJinv(0);
- // Reset matrices
- bUp = Eigen::VectorXd::Zero(nRow);
- bLw = Eigen::VectorXd::Zero(nRow);
-
- // Volume contribution
- // V_inf & dN' (stabilization)
- Eigen::MatrixXd dN(nRow, nDim);
- for (size_t i = 0; i < nRow; ++i)
- for (size_t j = 0; j < nDim; ++j)
- dN(i, j) = volUpDff(j, vsMap.at(wkN[i].first));
- Eigen::VectorXd JdNV(nRow);
- JdNV = dN * volUpJ.transpose() * Vi;
-
- // V2
- double grad2Up = volUpE->computeGrad2(u, 0);
- double grad2Lw = volLwE->computeGrad2(u, 0);
-
- // Build
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // N
- Eigen::VectorXd const &ff = surCache.getSf(k);
- Eigen::VectorXd N(nRow);
- for (size_t i = 0; i < nRow; ++i)
- N(i) = 0.5 * sqrt(surMem.getVol()) * JdNV(i) + ff(wkN[i].first);
-
- // Compute flux vector
- bUp += N * grad2Up * surGauss.getW(k) * surMem.getDetJ(k);
- bLw += N * grad2Lw * surGauss.getW(k) * surMem.getDetJ(k);
- }
-}
-
-/**
- * @brief Compute gradient of residual with respect to angle of attack
- *
- * s_i = int V*V * dV_inf*dN_i dS
- */
-void WakeElement::buildGradientResidualAoA(Eigen::VectorXd &bUp, Eigen::VectorXd &bLw, std::vector<double> const &u, Eigen::VectorXd const &dVi)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surUpE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- Eigen::MatrixXd const &volUpDff = volUpE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &surMem = surUpE->getVMem();
- Eigen::MatrixXd const &volUpJ = volUpE->getVMem().getJinv(0);
- // Reset matrices
- bUp = Eigen::VectorXd::Zero(nRow);
- bLw = Eigen::VectorXd::Zero(nRow);
-
- // Volume contribution
- // V_inf & dN' (stabilization)
- Eigen::MatrixXd dN(nRow, nDim);
- for (size_t i = 0; i < nRow; ++i)
- for (size_t j = 0; j < nDim; ++j)
- dN(i, j) = volUpDff(j, vsMap.at(wkN[i].first));
- Eigen::VectorXd JdNV(nRow);
- JdNV = dN * volUpJ.transpose() * dVi;
- // V2
- double grad2Up = volUpE->computeGrad2(u, 0);
- double grad2Lw = volLwE->computeGrad2(u, 0);
-
- // Build
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // N
- Eigen::VectorXd N(nRow);
- for (size_t i = 0; i < nRow; ++i)
- N(i) = 0.5 * sqrt(surMem.getVol()) * JdNV(i);
-
- // Compute flux vector
- bUp += N * grad2Up * surGauss.getW(k) * surMem.getDetJ(k);
- bLw += N * grad2Lw * surGauss.getW(k) * surMem.getDetJ(k);
- }
-}
-
-/**
- * @brief Compute gradient of residual with respect to angle of attack
- *
- * K_ij = int d(2*V*dN_j * V_inf*dN_i) dS
- * = sum 2*V*dgrad_phi * V_inf*gradN_i * detJ
- * + sum grad_phi^2 * dgradN_i detJ
- * + sum grad_phi^2 * V_inf*gradN_i ddetJ
- */
-void WakeElement::buildGradientResidualMesh(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, Eigen::MatrixXd &Ks, std::vector<double> const &u, Eigen::VectorXd const &vi)
-{
- // Get shape functions and Gauss points
- Cache &surCache = surUpE->getVCache();
- Gauss &surGauss = surCache.getVGauss();
- Eigen::MatrixXd const &dSfUp = volUpE->getVCache().getDsf(0);
- // Get Jacobian
- Mem &surMem = surUpE->getVMem();
- Mem &volUpMem = volUpE->getVMem();
- Mem &volLwMem = volLwE->getVMem();
- Eigen::MatrixXd const &iJUp = volUpMem.getJinv(0);
- Eigen::MatrixXd const &iJLw = volLwMem.getJinv(0);
- std::vector<Eigen::MatrixXd> const &dJUp = volUpMem.getGradJ(0);
- std::vector<Eigen::MatrixXd> const &dJLw = volLwMem.getGradJ(0);
- // Get surface gradient
- std::vector<double> const &dSurf = surMem.getGradVol();
- // Reset matrices
- Kup = Eigen::MatrixXd::Zero(nRow, nDim * volUpE->nodes.size());
- Klw = Eigen::MatrixXd::Zero(nRow, nDim * volLwE->nodes.size());
- Ks = Eigen::MatrixXd::Zero(nRow, nDim * surUpE->nodes.size());
-
- // Stabilization term
- double sqSurf = sqrt(surMem.getVol());
- Eigen::MatrixXd gradN(nRow, nDim);
- for (size_t i = 0; i < nRow; ++i)
- for (size_t m = 0; m < nDim; ++m)
- gradN(i, m) = dSfUp(m, vsMap.at(wkN[i].first));
- Eigen::VectorXd JgradNV = 0.5 * sqSurf * gradN * iJUp.transpose() * vi;
- // Gradient of stabilization term
- std::vector<Eigen::VectorXd> dGradSfV(nDim * volUpE->nodes.size());
- std::vector<Eigen::VectorXd> dGradSfS(nDim * surUpE->nodes.size());
- for (size_t i = 0; i < volUpE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- dGradSfV[idx] = 0.5 * sqSurf * gradN * (-iJUp * dJUp[idx] * iJUp).transpose() * vi; // S * dgradN
- }
- }
- for (size_t i = 0; i < surUpE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * gradN * iJUp.transpose() * vi; // dS * gradN
- }
- }
- // Gradients of potential
- Eigen::VectorXd dPhiUp = volUpE->computeGrad(u, 0);
- Eigen::VectorXd dPhiLw = volLwE->computeGrad(u, 0);
-
- for (size_t k = 0; k < surGauss.getN(); ++k)
- {
- // weight and Jacobian determinant
- double w = surGauss.getW(k);
- double detJ = surMem.getDetJ(k);
- std::vector<double> dDetJ;
- surUpE->buildGradientJ(k, dDetJ);
- // shape functions and stabilization
- Eigen::VectorXd const &sf = surCache.getSf(k);
- Eigen::VectorXd sfp(nRow);
- for (size_t i = 0; i < nRow; ++i)
- sfp(i) = sf(wkN[i].first) + JgradNV(i);
-
- // Build gradient wrt to upper volume nodes
- for (size_t i = 0; i < volUpE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- // derivative of velocity jump
- Kup.col(idx) += w * sfp * 2 * dPhiUp.transpose() * (-iJUp * dJUp[idx] * dPhiUp) * detJ; // (N + gradN') * dgrad_phi^2 * detJ
- // derivative of stabilization term (volume)
- Kup.col(idx) += w * dGradSfV[idx] * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * detJ; // dgradN' * grad_phi^2 * detJ
- }
- }
- // Build gradient wrt to lower volume nodes
- for (size_t i = 0; i < volLwE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- // derivative of velocity jump
- Klw.col(idx) += w * sfp * 2 * dPhiLw.transpose() * (-iJLw * dJLw[idx] * dPhiLw) * detJ; // (N + gradN') * dgrad_phi^2 * detJ
- }
- }
- // Build gradient wrt to (upper) surface nodes
- for (size_t i = 0; i < surUpE->nodes.size(); ++i)
- {
- for (size_t m = 0; m < nDim; ++m)
- {
- size_t idx = i * nDim + m;
- // derivative of stabilization term (surface)
- Ks.col(idx) += w * dGradSfS[idx] * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * detJ; // dgradN' * grad_phi^2 * detJ
- // derivative of jacobian determinant
- Ks.col(idx) += w * sfp * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * dDetJ[idx]; // (N + gradN') * grad_phi^2 * ddetJ
- }
- }
- }
-}
-
void WakeElement::write(std::ostream &out) const
{
out << "WakeElement #" << no
diff --git a/flow/src/wWakeElement.h b/flow/src/wWakeElement.h
index 28ba1347..cf531650 100644
--- a/flow/src/wWakeElement.h
+++ b/flow/src/wWakeElement.h
@@ -47,26 +47,21 @@ public:
tbox::Element *volUpE; ///< Connected volume elements "+" side
tbox::Element *volLwE; ///< Connected volume elements "-" side
- size_t nRow; ///< number of rows in stiffness matrix
- size_t nColUp; ///< number of columns in in stiffness matrix "+" side
- size_t nColLw; ///< number of columns in in stiffness matrix "-" side
+ size_t nRow; ///< number of rows in stiffness matrix
+ size_t nColUp; ///< number of columns in in stiffness matrix "+" side
+ size_t nColLw; ///< number of columns in in stiffness matrix "-" side
+ size_t nDim; ///< dimension of volume element
+
std::vector<std::pair<size_t, tbox::Node *>> wkN; ///< map between local index and contributing wake global nodes
+ std::map<size_t, size_t> vsMap; ///< map between surface node and volume node indices
WakeElement(size_t _no, tbox::Element *&_surUpE, tbox::Element *&_surLwE, tbox::Element *&_volUpE, tbox::Element *&_volLwE, std::vector<std::pair<size_t, tbox::Node *>> &_wkN);
#ifndef SWIG
- void buildK(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, std::vector<double> const &u, Eigen::VectorXd const &Vi);
- void buildNK(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, std::vector<double> const &u, Eigen::VectorXd const &Vi);
- void buildNs(Eigen::VectorXd &bUp, Eigen::VectorXd &bLw, std::vector<double> const &u, Eigen::VectorXd const &Vi);
- void buildGradientResidualAoA(Eigen::VectorXd &bUp, Eigen::VectorXd &bLw, std::vector<double> const &u, Eigen::VectorXd const &dVi);
- void buildGradientResidualMesh(Eigen::MatrixXd &Kup, Eigen::MatrixXd &Klw, Eigen::MatrixXd &Ks, std::vector<double> const &u, Eigen::VectorXd const &vi);
virtual void write(std::ostream &out) const override;
#endif
private:
- size_t nDim; ///< dimension of volume element
- std::map<size_t, size_t> vsMap; ///< map between surface node and volume node indices
-
void mapNodes();
};
diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp
new file mode 100644
index 00000000..376db08f
--- /dev/null
+++ b/flow/src/wWakeResidual.cpp
@@ -0,0 +1,251 @@
+/*
+ * 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 "wWakeResidual.h"
+#include "wWakeElement.h"
+
+#include "wElement.h"
+#include "wCache.h"
+#include "wGauss.h"
+#include "wMem.h"
+
+using namespace tbox;
+using namespace flow;
+
+/**
+ * @brief Build the residual matrix for a fixed-point iteration, on one wake element
+ *
+ * A_u = \int (psi + grad_psi*v_inf) * v_u*grad_phi_u dS
+ * A_l = \int (psi + grad_psi*v_inf) * v_l*grad_phi_l dS
+ * NB: A_u and A_l will be assembled on same row DOF, but on different columns
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+{
+ // Get shape functions and Gauss points
+ Cache &sCache = we.surUpE->getVCache();
+ Gauss &sGauss = sCache.getVGauss();
+ Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
+ Eigen::MatrixXd const &vLwDsf = we.volLwE->getVCache().getDsf(0);
+ // Get Jacobian
+ Mem &sMem = we.surUpE->getVMem();
+ Eigen::MatrixXd const &vUpJ = we.volUpE->getVMem().getJinv(0);
+ Eigen::MatrixXd const &vLwJ = we.volLwE->getVMem().getJinv(0);
+ // Get size
+ size_t nRow = we.nRow;
+ size_t nDim = we.nDim;
+
+ // Stabilization
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t j = 0; j < nDim; ++j)
+ dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
+ Eigen::VectorXd dSfp(nRow);
+ dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * vUpJ.transpose() * vi;
+ // Velocity
+ Eigen::RowVectorXd Aup = we.volUpE->computeGrad(phi, 0).transpose() * vUpJ * vUpDsf;
+ Eigen::RowVectorXd Alw = we.volLwE->computeGrad(phi, 0).transpose() * vLwJ * vLwDsf;
+
+ // Build
+ Eigen::MatrixXd Au = Eigen::MatrixXd::Zero(nRow, we.nColUp);
+ Eigen::MatrixXd Al = Eigen::MatrixXd::Zero(nRow, we.nColLw);
+ for (size_t k = 0; k < sGauss.getN(); ++k)
+ {
+ // stabilized shape functions
+ Eigen::VectorXd const &sf = sCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = dSfp(i) + sf(we.wkN[i].first);
+ // wake contribution
+ Au += sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
+ Al += sfp * Alw * sGauss.getW(k) * sMem.getDetJ(k);
+ }
+ return std::make_tuple(Au, Al);
+}
+
+/**
+ * @brief Build the residual vector, on one wake element
+ *
+ * b_u = \int (psi + grad_psi*v_inf) * (grad_phi_u)^2 dS
+ * b_l = \int (psi + grad_psi*v_inf) * (grad_phi_l)^2 dS
+ */
+std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+{
+ // Get shape functions and Gauss points
+ Cache &sCache = we.surUpE->getVCache();
+ Gauss &sGauss = sCache.getVGauss();
+ Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
+ // Get Jacobian
+ Mem &sMem = we.surUpE->getVMem();
+ // Get size
+ size_t nRow = we.nRow;
+ size_t nDim = we.nDim;
+
+ // Stabilization
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t j = 0; j < nDim; ++j)
+ dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
+ Eigen::VectorXd dSfp(nRow);
+ dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * we.volUpE->getVMem().getJinv(0).transpose() * vi;
+ // Velocity squared
+ double dPhi2u = we.volUpE->computeGrad2(phi, 0);
+ double dPhi2l = we.volLwE->computeGrad2(phi, 0);
+
+ // Build
+ Eigen::VectorXd bu = Eigen::VectorXd::Zero(nRow);
+ Eigen::VectorXd bl = Eigen::VectorXd::Zero(nRow);
+ for (size_t k = 0; k < sGauss.getN(); ++k)
+ {
+ // stabilized shape functions
+ Eigen::VectorXd const &sf = sCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = dSfp(i) + sf(we.wkN[i].first);
+ // wake contribution
+ bu += sfp * dPhi2u * sGauss.getW(k) * sMem.getDetJ(k);
+ bl += sfp * dPhi2l * sGauss.getW(k) * sMem.getDetJ(k);
+ }
+ return std::make_tuple(bu, bl);
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one wake element
+ *
+ * A_u = \int (psi + grad_psi*v_inf) * 2*(grad_phi_u)*dgrad_phi_u dS
+ * A_l = \int (psi + grad_psi*v_inf) * 2*(grad_phi_l)*dgrad_phi_l dS
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+{
+ // Build
+ Eigen::MatrixXd Au, Al;
+ std::tie(Au, Al) = WakeResidual::buildFixed(we, phi, vi);
+ return std::make_tuple(2* Au, 2* Al);
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the nodes, on one wake element
+ *
+ * A = d( \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS )
+ * = \int dgrad_psi*v_inf * (grad_phi)^2 dS
+ * + \int (psi + grad_psi*v_inf) * d(grad_phi)^2 dS
+ * + \int (psi + grad_psi*v_inf) * (grad_phi)^2 ddS
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+{
+ // Get shape functions and Gauss points
+ Cache &sCache = we.surUpE->getVCache();
+ Gauss &surGauss = sCache.getVGauss();
+ Eigen::MatrixXd const &dSfUp = we.volUpE->getVCache().getDsf(0);
+ // Get Jacobian
+ Mem &sMem = we.surUpE->getVMem();
+ Mem &vUpMem = we.volUpE->getVMem();
+ Mem &vLwMem = we.volLwE->getVMem();
+ Eigen::MatrixXd const &iJUp = vUpMem.getJinv(0);
+ Eigen::MatrixXd const &iJLw = vLwMem.getJinv(0);
+ std::vector<Eigen::MatrixXd> const &dJUp = vUpMem.getGradJ(0);
+ std::vector<Eigen::MatrixXd> const &dJLw = vLwMem.getGradJ(0);
+ // Get surface gradient
+ std::vector<double> const &dSurf = sMem.getGradVol();
+ // Get size
+ size_t nRow = we.nRow;
+ size_t nDim = we.nDim;
+ size_t nColUp = we.nColUp;
+ size_t nColLw = we.nColLw;
+ size_t nSur = we.surUpE->nodes.size();
+
+ // Stabilization term
+ double sqSurf = sqrt(sMem.getVol());
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t m = 0; m < nDim; ++m)
+ dSf(i, m) = dSfUp(m, we.vsMap.at(we.wkN[i].first));
+ Eigen::VectorXd gradSfp = 0.5 * sqSurf * dSf * iJUp.transpose() * vi;
+ // Gradient of stabilization term (volume and surface)
+ std::vector<Eigen::VectorXd> dGradSfV(nDim * nColUp);
+ for (size_t i = 0; i < nColUp; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJUp * dJUp[idx] * iJUp).transpose() * vi; // S * dgrad_psi
+ }
+ }
+ std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur);
+ for (size_t i = 0; i < nSur; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJUp.transpose() * vi; // dS * grad_psi
+ }
+ }
+ // Gradients of potential
+ Eigen::VectorXd dPhiUp = we.volUpE->computeGrad(phi, 0);
+ Eigen::VectorXd dPhiLw = we.volLwE->computeGrad(phi, 0);
+
+ // Build
+ Eigen::MatrixXd Au = Eigen::MatrixXd::Zero(nRow, nDim * nColUp);
+ Eigen::MatrixXd Al = Eigen::MatrixXd::Zero(nRow, nDim * nColLw);
+ Eigen::MatrixXd As = Eigen::MatrixXd::Zero(nRow, nDim * nSur);
+ for (size_t k = 0; k < surGauss.getN(); ++k)
+ {
+ // weight and Jacobian determinant
+ double w = surGauss.getW(k);
+ double detJ = sMem.getDetJ(k);
+ std::vector<double> dDetJ = sMem.getGradDetJ(k);
+ // shape functions and stabilization
+ Eigen::VectorXd const &sf = sCache.getSf(k);
+ Eigen::VectorXd sfp(nRow);
+ for (size_t i = 0; i < nRow; ++i)
+ sfp(i) = sf(we.wkN[i].first) + gradSfp(i);
+
+ // Build gradient wrt to upper volume nodes
+ for (size_t i = 0; i < nColUp; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ // gradient of velocity jump
+ Au.col(idx) += w * sfp * 2 * dPhiUp.transpose() * (-iJUp * dJUp[idx] * dPhiUp) * detJ; // (psi + grad_psi) * dgrad_phi^2 * detJ
+ // gradient of stabilization term (volume)
+ Au.col(idx) += w * dGradSfV[idx] * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * detJ; // dgrad_psi' * grad_phi^2 * detJ
+ }
+ }
+ // Build gradient wrt to lower volume nodes
+ for (size_t i = 0; i < nColLw; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ // gradient of velocity jump
+ Al.col(idx) += w * sfp * 2 * dPhiLw.transpose() * (-iJLw * dJLw[idx] * dPhiLw) * detJ; // (psi + grad_psi) * dgrad_phi^2 * detJ
+ }
+ }
+ // Build gradient wrt to (upper) surface nodes
+ for (size_t i = 0; i < nSur; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ // gradient of stabilization term (surface)
+ As.col(idx) += w * dGradSfS[idx] * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * detJ; // dgrad_psi * grad_phi^2 * detJ
+ // gradient of jacobian determinant
+ As.col(idx) += w * sfp * (dPhiUp.dot(dPhiUp) - dPhiLw.dot(dPhiLw)) * dDetJ[idx]; // (N + grad_psi) * grad_phi^2 * ddetJ
+ }
+ }
+ }
+ return std::make_tuple(Au, Al, As);
+}
diff --git a/flow/src/wWakeResidual.h b/flow/src/wWakeResidual.h
new file mode 100644
index 00000000..fc536849
--- /dev/null
+++ b/flow/src/wWakeResidual.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.
+ */
+
+#ifndef WWAKERESIDUAL_H
+#define WWAKERESIDUAL_H
+
+#include "flow.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Formulation of wake boundary residuals
+ */
+class FLOW_API WakeResidual
+{
+public:
+ // Picard
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildFixed(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi);
+ // Newton
+ static std::tuple<Eigen::VectorXd, Eigen::VectorXd> build(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi);
+ // Adjoint
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi);
+};
+
+} // namespace flow
+#endif //WWAKERESIDUAL_H
diff --git a/flow/tests/lift.py b/flow/tests/lift.py
index 5a9f4959..a5757ec9 100644
--- a/flow/tests/lift.py
+++ b/flow/tests/lift.py
@@ -89,13 +89,16 @@ def main():
if M_inf == 0 and alpha == 3*math.pi/180:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.1, 1.5e-1))
tests.add(CTest('Cl', solver.Cl, 0.36, 5e-2))
+ tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
elif M_inf == 0.6 and alpha == 2*math.pi/180:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.03, 1e-1))
tests.add(CTest('Cl', solver.Cl, 0.315, 5e-2))
+ tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
elif M_inf == 0.7 and alpha == 2*math.pi/180:
tests.add(CTest('iteration count', solver.nIt, 12, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.28, 5e-2))
tests.add(CTest('Cl', solver.Cl, 0.388, 5e-2))
+ tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
diff --git a/flow/tests/lift3.py b/flow/tests/lift3.py
index cc0bd631..da872d0f 100644
--- a/flow/tests/lift3.py
+++ b/flow/tests/lift3.py
@@ -89,6 +89,8 @@ def main():
tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
+ tests.add(CTest('CS', solver.Cs, 0.012, 1e-2))
+ tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
diff --git a/flow/tests/meshDef.py b/flow/tests/meshDef.py
index 411b0c6b..d258e5a6 100644
--- a/flow/tests/meshDef.py
+++ b/flow/tests/meshDef.py
@@ -64,7 +64,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, M_crit, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, M_crit, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te')
# set the refrence problem and add medium, initial/boundary conditions
pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, M_crit, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te')
tms['pre'].stop()
@@ -121,7 +121,6 @@ def main():
print(' case M alpha Cl Cd Cm')
print(' true {0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, solver.Cm))
print(' ref {0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alfa*180/np.pi, solver_ref.Cl, solver_ref.Cd, solver_ref.Cm))
- print('Cm will differ unless dx and dy are set to 0...')
# display timers
tms["total"].stop()
@@ -140,6 +139,7 @@ def main():
tests.add(CTest('solver_ref: iteration count', solver_ref.nIt, 5, 1, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), min(Cp_ref[:,3]), 5e-2))
tests.add(CTest('Cl', solver.Cl, solver_ref.Cl, 5e-2))
+ tests.add(CTest('Cm', solver.Cm, solver_ref.Cm, 5e-2))
tests.run()
# eof
diff --git a/flow/tests/meshDef3.py b/flow/tests/meshDef3.py
index 1b2d9502..ab3fd22a 100644
--- a/flow/tests/meshDef3.py
+++ b/flow/tests/meshDef3.py
@@ -124,6 +124,8 @@ def main():
tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.135, 5e-1))
tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
+ tests.add(CTest('CS', solver.Cs, -0.011, 1e-2))
+ tests.add(CTest('CM', solver.Cm, 0.0061, 1e-2))
tests.run()
# eof
diff --git a/flow/tests/nonlift.py b/flow/tests/nonlift.py
index 450008c8..ef6c22cc 100644
--- a/flow/tests/nonlift.py
+++ b/flow/tests/nonlift.py
@@ -96,7 +96,8 @@ def main():
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
else:
raise Exception('Test not defined for this flow')
- tests.add(CTest('Cl', solver.Cl, 0., 5e-2))
+ tests.add(CTest('Cl', solver.Cl, 0., 1e-2))
+ tests.add(CTest('Cm', solver.Cm, 0., 1e-2))
tests.run()
# eof
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index a1b3d561..303f6efa 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -205,10 +205,6 @@ Eigen::MatrixXd Element::buildK(std::vector<double> const &u, Fct2 const &fct, b
{
throw std::runtime_error("Element::buildK not implemented\n");
}
-Eigen::MatrixXd Element::buildDK(std::vector<double> const &u, Fct0 const &fct)
-{
- throw std::runtime_error("Element::buildNK not implemented\n");
-}
//---------------------
Eigen::MatrixXd Element::buildK_mechanical(Eigen::MatrixXd const &H)
{
@@ -265,19 +261,15 @@ Eigen::VectorXd Element::builds(std::vector<double> const &u, Fct1 const &f)
{
throw std::runtime_error("Element::builds not implemented\n");
}
-Eigen::VectorXd Element::buildDs(std::vector<double> const &u, Fct0 const &f)
-{
- throw std::runtime_error("Element::buildNs not implemented\n");
-}
Eigen::MatrixXd Element::build(MATTYPE type)
{
throw std::runtime_error("Element::build not implemented\n");
}
-Eigen::VectorXd Element::computeGrad(std::vector<double> const &u, size_t k)
+Eigen::VectorXd Element::computeGrad(std::vector<double> const &u, size_t k) const
{
throw std::runtime_error("Element::computeGrad not implemented\n");
}
-double Element::computeGrad2(std::vector<double> const &u, size_t k)
+double Element::computeGrad2(std::vector<double> const &u, size_t k) const
{
throw std::runtime_error("Element::computeGrad2 not implemented\n");
}
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 48d60c92..068f1bec 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -94,7 +94,6 @@ public:
virtual void write(std::ostream &out) const override;
virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct);
virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct2 const &fct, bool fake);
- virtual Eigen::MatrixXd buildDK(std::vector<double> const &u, Fct0 const &fct);
//---------------------
virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H);
virtual Eigen::MatrixXd buildS_thermomechanical(Eigen::MatrixXd const &D);
@@ -109,10 +108,9 @@ public:
virtual Eigen::MatrixXd buildS();
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f);
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f);
- virtual Eigen::VectorXd buildDs(std::vector<double> const &u, Fct0 const &f);
virtual Eigen::MatrixXd build(MATTYPE type);
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k);
- virtual double computeGrad2(std::vector<double> const &u, size_t k);
+ virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const;
+ virtual double computeGrad2(std::vector<double> const &u, size_t k) const;
virtual Eigen::VectorXd computeqV(std::vector<double> const &u, Fct2 const &fct);
virtual Eigen::VectorXd computeqint(std::vector<double> const &u, Fct2 const &fct);
virtual double computeV(std::vector<double> const &u, Fct0 const &fct);
diff --git a/tbox/src/wFct0.cpp b/tbox/src/wFct0.cpp
index 4b3baad3..b0a7b0dc 100644
--- a/tbox/src/wFct0.cpp
+++ b/tbox/src/wFct0.cpp
@@ -66,7 +66,7 @@ void PwLf::write(std::ostream &out) const
//-------------------------------------------------------------------------
-double Fct0::evalD(Element &e, std::vector<double> const &u, size_t k) const
+double Fct0::evalD(Element const &e, std::vector<double> const &u, size_t k) const
{
throw std::runtime_error("Fct0::evalD not implemented\n");
}
@@ -81,7 +81,7 @@ void Fct0C::write(std::ostream &out) const
//-------------------------------------------------------------------------
double
-Fct0XYZ::eval(Element &e, std::vector<double> const &u, size_t k) const
+Fct0XYZ::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
Cache &cache = e.getVCache();
Eigen::VectorXd const &ff = cache.getSf(k);
@@ -102,7 +102,7 @@ void Fct0XYZ::write(std::ostream &out) const
//-------------------------------------------------------------------------
double
-Fct0U::eval(Element &e, std::vector<double> const &u, size_t k) const
+Fct0U::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
Cache &cache = e.getVCache();
Eigen::VectorXd const &ff = cache.getSf(k);
diff --git a/tbox/src/wFct0.h b/tbox/src/wFct0.h
index a46d59cb..81da6cd4 100644
--- a/tbox/src/wFct0.h
+++ b/tbox/src/wFct0.h
@@ -59,10 +59,10 @@ public:
}
#ifndef SWIG
- virtual double eval(Element &e,
+ virtual double eval(Element const &e,
std::vector<double> const &u,
size_t k) const = 0;
- virtual double evalD(Element &e,
+ virtual double evalD(Element const &e,
std::vector<double> const &u,
size_t k) const;
#endif
@@ -81,13 +81,13 @@ public:
void update(double _v) { v = _v; }
#ifndef SWIG
- virtual double eval(Element &e,
+ virtual double eval(Element const &e,
std::vector<double> const &u,
size_t k) const override
{
return v;
}
- virtual double evalD(Element &e,
+ virtual double evalD(Element const &e,
std::vector<double> const &u,
size_t k) const override
{
@@ -107,7 +107,7 @@ public:
virtual double eval(Eigen::Vector3d const &pos) const = 0;
#ifndef SWIG
- virtual double eval(Element &e,
+ virtual double eval(Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
@@ -124,7 +124,7 @@ public:
virtual double eval(double u) const = 0;
#ifndef SWIG
- virtual double eval(Element &e,
+ virtual double eval(Element const &e,
std::vector<double> const &u,
size_t k) const override;
virtual void write(std::ostream &out) const override;
diff --git a/tbox/src/wFct1.h b/tbox/src/wFct1.h
index 0edd8584..439747f2 100644
--- a/tbox/src/wFct1.h
+++ b/tbox/src/wFct1.h
@@ -44,9 +44,12 @@ public:
}
#ifndef SWIG
- virtual Eigen::Vector3d eval(Element &e,
+ virtual Eigen::Vector3d eval(Element const &e,
std::vector<double> const &u,
size_t k) const = 0;
+ virtual Eigen::Vector3d evalD(tbox::Element const &e,
+ std::vector<double> const &u,
+ size_t k) const = 0;
#endif
};
@@ -67,12 +70,18 @@ public:
}
#ifndef SWIG
- virtual Eigen::Vector3d eval(Element &e,
+ virtual Eigen::Vector3d eval(Element const &e,
std::vector<double> const &u,
size_t k) const override
{
return v;
}
+ virtual Eigen::Vector3d evalD(tbox::Element const &e,
+ std::vector<double> const &u,
+ size_t k) const override
+ {
+ return Eigen::Vector3d::Zero();
+ }
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/tbox/src/wFct2.cpp b/tbox/src/wFct2.cpp
index b510d114..4b83977a 100644
--- a/tbox/src/wFct2.cpp
+++ b/tbox/src/wFct2.cpp
@@ -35,7 +35,7 @@ Fct2C::Fct2C(double _v11, double _v22, double _v12) : Fct2()
throw std::runtime_error("det(Fct2C)<0!");
}
-void Fct2C::eval(Element &e,
+void Fct2C::eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const
{
@@ -49,7 +49,7 @@ void Fct2C::write(std::ostream &out) const
//-------------------------------------------------------------------------
-void Fct2XYZ::eval(Element &e,
+void Fct2XYZ::eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const
{
@@ -71,7 +71,7 @@ void Fct2XYZ::write(std::ostream &out) const
//-------------------------------------------------------------------------
-void Fct2U::eval(Element &e,
+void Fct2U::eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const
{
@@ -93,7 +93,7 @@ void Fct2U::write(std::ostream &out) const
//-------------------------------------------------------------------------
-void Fct2UdU::eval(Element &e,
+void Fct2UdU::eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const
{
diff --git a/tbox/src/wFct2.h b/tbox/src/wFct2.h
index 002a2188..4b902552 100644
--- a/tbox/src/wFct2.h
+++ b/tbox/src/wFct2.h
@@ -45,7 +45,7 @@ public:
virtual void evalall() {}
#ifndef SWIG
- virtual void eval(Element &e,
+ virtual void eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const = 0;
#endif
@@ -62,7 +62,7 @@ public:
Fct2C(double _v11, double _v22, double _v12);
#ifndef SWIG
- virtual void eval(Element &e,
+ virtual void eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const override;
virtual void write(std::ostream &out) const override;
@@ -78,7 +78,7 @@ public:
virtual void eval(Eigen::Vector3d const &pos, Eigen::MatrixXd &v, bool fake) const = 0;
#ifndef SWIG
- virtual void eval(Element &e,
+ virtual void eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const override;
virtual void write(std::ostream &out) const override;
@@ -94,7 +94,7 @@ public:
virtual void eval(double u, Eigen::MatrixXd &v, bool fake) const = 0;
#ifndef SWIG
- virtual void eval(Element &e,
+ virtual void eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const override;
virtual void write(std::ostream &out) const override;
@@ -106,10 +106,10 @@ public:
class TBOX_API Fct2UdU : public Fct2
{
public:
- virtual void eval(Element &e, size_t k, double u, Eigen::VectorXd const &gradu,
+ virtual void eval(Element const &e, size_t k, double u, Eigen::VectorXd const &gradu,
Eigen::MatrixXd &v, bool fake) const = 0;
#ifndef SWIG
- virtual void eval(Element &e,
+ virtual void eval(Element const &e,
std::vector<double> const &u,
size_t k, Eigen::MatrixXd &out, bool fake) const override;
virtual void write(std::ostream &out) const override;
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 986904ef..9e34ca22 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -203,9 +203,9 @@ Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct1 const &f)
*
* grad = (Delta_U / Delta_L)^2
*/
-double Line2::computeGrad2(std::vector<double> const &u, size_t k)
+double Line2::computeGrad2(std::vector<double> const &u, size_t k) const
{
- double dUdL = (u[nodes[1]->row] - u[nodes[0]->row]) / computeV(); // V = length of the element
+ double dUdL = (u[nodes[1]->row] - u[nodes[0]->row]) / getMem().getVol(); // V = length of the element
return dUdL * dUdL;
}
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 400d87da..4e033f3e 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -43,7 +43,7 @@ class TBOX_API Line2 : public Element
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f) override;
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f) override;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) override;
+ virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
virtual double computeV() override;
diff --git a/tbox/src/wMem.inl b/tbox/src/wMem.inl
index 17aa9331..217af82a 100644
--- a/tbox/src/wMem.inl
+++ b/tbox/src/wMem.inl
@@ -97,7 +97,7 @@ inline std::vector<double> const &Mem::getGradDetJ(size_t k)
}
/**
- * @brief Return the gradient of determinant of the Jacobian with respect to node positions at Gauss point k
+ * @brief Return the gradient of the Jacobian matrix with respect to node positions at Gauss point k
*/
inline std::vector<Eigen::MatrixXd> const &Mem::getGradJ(size_t k)
{
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index fe2e1fca..e488be58 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -165,50 +165,6 @@ Eigen::MatrixXd Tetra4::buildK(std::vector<double> const &u, Fct0 const &fct)
return K;
}
-/**
- * @brief Compute (first part of) derivative stiffness matrix for scalar field (Newton method)
- *
- * K_ij(4,4) = int (df du*dN_j du*dN_i) dV
- */
-Eigen::MatrixXd Tetra4::buildDK(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd K = Eigen::Matrix4d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse, shape functions and gradient
- Eigen::Matrix3d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::VectorXd gradu = computeGrad(u, k);
-
- // Elementary stiffness matrix
- K += (gradu.transpose() * J * dff).transpose() * gradu.transpose() * J * dff * fct.evalD(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
-/**
- * @brief Compute derivative of flux vector for a scalar field (Newton method)
- *
- * s_i(4) = int f*du*dN_i dV
- */
-Eigen::VectorXd Tetra4::buildDs(std::vector<double> const &u, Fct0 const &f)
-{
- CacheTetra4 &cache = getCache(order);
- GaussTetra4 &gauss = cache.gauss;
- MemTetra4 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector4d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += computeGrad(u, k).transpose() * mem.getJinv(k) * cache.getDsf(k) * f.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k); // elementary flux vector
- return s;
-}
-
/**
* @brief Compute linear elasticity stiffness matrix
*
@@ -337,7 +293,7 @@ void Tetra4::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eig
/**
* @brief Compute gradient of a scalar field at Gauss point k
*/
-Eigen::VectorXd Tetra4::computeGrad(std::vector<double> const &u, size_t k)
+Eigen::VectorXd Tetra4::computeGrad(std::vector<double> const &u, size_t k) const
{
// Solution vector
Eigen::Vector4d ue;
@@ -350,7 +306,7 @@ Eigen::VectorXd Tetra4::computeGrad(std::vector<double> const &u, size_t k)
/**
* @brief Compute square of gradient of a scalar field at Gauss point k
*/
-double Tetra4::computeGrad2(std::vector<double> const &u, size_t k)
+double Tetra4::computeGrad2(std::vector<double> const &u, size_t k) const
{
return computeGrad(u, k).squaredNorm();
}
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index 76108f23..f524f9c9 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -42,16 +42,14 @@ class TBOX_API Tetra4 : public Element
virtual Eigen::MatrixXd buildM() override;
virtual Eigen::MatrixXd buildK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::MatrixXd buildDK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::VectorXd buildDs(std::vector<double> const &u, Fct0 const &f) override;
virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H) override;
virtual Eigen::MatrixXd buildS_thermomechanical(Eigen::MatrixXd const &D) override;
virtual Eigen::MatrixXd buildL_thermal(Eigen::MatrixXd const &C) override;
virtual void strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u) override;
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) override;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) override;
+ virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const override;
+ virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
virtual double computeV() override;
#endif
diff --git a/tbox/src/wTri3.cpp b/tbox/src/wTri3.cpp
index b5a107db..2be7fa46 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -318,32 +318,6 @@ Eigen::MatrixXd Tri3::buildK_mechanical(Eigen::MatrixXd const &H)
return K;
}
-/**
- * @brief Compute (first part of) derivative stiffness matrix for scalar field (Newton method)
- *
- * K_ij(3,3) = int (df du*dN_j du*dN_i) dV
- */
-Eigen::MatrixXd Tri3::buildDK(std::vector<double> const &u, Fct0 const &fct)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::MatrixXd K = Eigen::Matrix3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- {
- // Jacobian inverse, shape functions and gradient
- Eigen::Matrix2d const &J = mem.getJinv(k);
- Eigen::MatrixXd const &dff = cache.getDsf(k);
- Eigen::VectorXd gradu = computeGrad(u, k);
-
- // Elementary stiffness matrix
- K += (gradu.transpose() * J * dff).transpose() * gradu.transpose() * J * dff * fct.evalD(*this, u, k) * gauss.getW(k) * mem.getDetJ(k);
- }
- return K;
-}
-
/**
* @brief Compute flux vector for a scalar field
*
@@ -380,24 +354,6 @@ Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct1 const &f)
return s;
}
-/**
- * @brief Compute derivative of flux vector for a scalar field (Newton method)
- *
- * s_i(3) = int f*du*dN_i dV
- */
-Eigen::VectorXd Tri3::buildDs(std::vector<double> const &u, Fct0 const &f)
-{
- CacheTri3 &cache = getCache(order);
- GaussTri3 &gauss = cache.gauss;
- MemTri3 &mem = getMem();
-
- // Gauss integration
- Eigen::VectorXd s = Eigen::Vector3d::Zero();
- for (size_t k = 0; k < gauss.getN(); ++k)
- s += computeGrad(u, k).transpose() * mem.getJinv(k) * cache.getDsf(k) * f.eval(*this, u, k) * gauss.getW(k) * mem.getDetJ(k); // elementary flux vector
- return s;
-}
-
/**
* @brief Compute mass matrix
*
@@ -425,7 +381,7 @@ Eigen::MatrixXd Tri3::buildS()
/**
* @brief Compute gradient of a scalar field at Gauss point k
*/
-Eigen::VectorXd Tri3::computeGrad(std::vector<double> const &u, size_t k)
+Eigen::VectorXd Tri3::computeGrad(std::vector<double> const &u, size_t k) const
{
// Solution vector
Eigen::Vector3d ue;
@@ -438,7 +394,7 @@ Eigen::VectorXd Tri3::computeGrad(std::vector<double> const &u, size_t k)
/**
* @brief Compute square of gradient of a scalar field at Gauss point k
*/
-double Tri3::computeGrad2(std::vector<double> const &u, size_t k)
+double Tri3::computeGrad2(std::vector<double> const &u, size_t k) const
{
return computeGrad(u, k).squaredNorm();
}
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index 7255a33c..bea66750 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -47,12 +47,10 @@ class TBOX_API Tri3 : public Element
virtual Eigen::MatrixXd buildK_mechanical(Eigen::MatrixXd const &H) override;
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct0 const &f) override;
virtual Eigen::VectorXd builds(std::vector<double> const &u, Fct1 const &f) override;
- virtual Eigen::MatrixXd buildDK(std::vector<double> const &u, Fct0 const &fct) override;
- virtual Eigen::VectorXd buildDs(std::vector<double> const &u, Fct0 const &f) override;
virtual Eigen::MatrixXd buildS() override;
- virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) override;
- virtual double computeGrad2(std::vector<double> const &u, size_t k) override;
+ virtual Eigen::VectorXd computeGrad(std::vector<double> const &u, size_t k) const override;
+ virtual double computeGrad2(std::vector<double> const &u, size_t k) const override;
virtual Eigen::VectorXd computeqV(std::vector<double> const &u, Fct2 const &fct) override;
virtual Eigen::VectorXd computeqint(std::vector<double> const &u, Fct2 const &fct) override;
virtual double computeV(std::vector<double> const &u, Fct0 const &fct) override;
--
GitLab