From 74883c918252105071e6effec757e6790600bdb2 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Mon, 28 Jun 2021 17:44:01 +0200
Subject: [PATCH] Add Observer. Remove functions from python side and add as
observer of Problem. Add loads directions as functions. Minor changes.
---
flow/_src/floww.h | 4 -
flow/_src/floww.i | 26 ------
flow/benchmark/onera.py | 3 +-
flow/cases/coyote.py | 1 -
flow/cases/n0012.py | 1 -
flow/cases/n0012_3.py | 1 -
flow/cases/n64A410.py | 1 -
flow/cases/rae2822.py | 1 -
flow/cases/wbht.py | 1 -
flow/cases/wht.py | 1 -
flow/default.py | 17 ++--
flow/scripts/config.py | 23 ++---
flow/scripts/polar.py | 9 +-
flow/scripts/trim.py | 9 +-
flow/src/flow.h | 13 ++-
flow/src/wAdjoint.cpp | 64 +++----------
flow/src/wAssign.cpp | 31 ++++---
flow/src/wAssign.h | 16 ++--
flow/src/wBlowing.cpp | 6 +-
flow/src/wBlowing.h | 4 +-
flow/src/wBlowingResidual.cpp | 4 +-
flow/src/wBlowingResidual.h | 2 +-
flow/src/wF0El.cpp | 69 ++++++--------
flow/src/wF0El.h | 136 +++++++++++++--------------
flow/src/wF0Ps.cpp | 25 +----
flow/src/wF0Ps.h | 41 ++------
flow/src/wF1Ct.cpp | 154 +++++++++++++++++++++++++++++++
flow/src/wF1Ct.h | 99 ++++++++++++++++++++
flow/src/wF1El.cpp | 14 +--
flow/src/wF1El.h | 36 +++++---
flow/src/wFreestream.cpp | 12 ++-
flow/src/wFreestream.h | 12 +--
flow/src/wFreestreamResidual.cpp | 4 +-
flow/src/wLoadFunctional.cpp | 12 +--
flow/src/wLoadFunctional.h | 6 +-
flow/src/wMedium.cpp | 37 +++++++-
flow/src/wMedium.h | 21 +++--
flow/src/wNewton.cpp | 10 +-
flow/src/wPicard.cpp | 5 +-
flow/src/wPotentialResidual.cpp | 37 ++++----
flow/src/wPotentialResidual.h | 1 +
flow/src/wProblem.cpp | 112 +++++-----------------
flow/src/wProblem.h | 37 ++++----
flow/src/wSolver.cpp | 24 ++---
flow/src/wWakeResidual.cpp | 19 ++--
flow/src/wWakeResidual.h | 8 +-
flow/tests/adjoint.py | 3 +-
flow/tests/adjoint3.py | 3 +-
flow/tests/bli.py | 3 +-
flow/tests/cylinder.py | 3 +-
flow/tests/cylinder2D5.py | 3 +-
flow/tests/cylinder3.py | 3 +-
flow/tests/lift.py | 3 +-
flow/tests/lift3.py | 3 +-
flow/tests/meshDef.py | 5 +-
flow/tests/meshDef3.py | 3 +-
flow/tests/nonlift.py | 3 +-
flow/validation/agard.py | 4 +-
flow/validation/onera.py | 3 +-
fwk/src/fwk.h | 1 +
fwk/src/wObserver.cpp | 25 +++++
fwk/src/wObserver.h | 39 ++++++++
heat/src/wSolver.cpp | 1 -
tbox/src/wFct0.cpp | 7 --
tbox/src/wFct0.h | 15 +--
tbox/src/wFct1.h | 11 +--
66 files changed, 712 insertions(+), 598 deletions(-)
create mode 100644 flow/src/wF1Ct.cpp
create mode 100644 flow/src/wF1Ct.h
create mode 100644 fwk/src/wObserver.cpp
create mode 100644 fwk/src/wObserver.h
diff --git a/flow/_src/floww.h b/flow/_src/floww.h
index b593315c..297a990a 100644
--- a/flow/_src/floww.h
+++ b/flow/_src/floww.h
@@ -18,11 +18,7 @@
#include "wAssign.h"
#include "wBlowing.h"
#include "wBoundary.h"
-#include "wFace.h"
#include "wFreestream.h"
-#include "wF0El.h"
-#include "wF0Ps.h"
-#include "wF1El.h"
#include "wKutta.h"
#include "wKuttaElement.h"
#include "wMedium.h"
diff --git a/flow/_src/floww.i b/flow/_src/floww.i
index 4f7387e4..3e969f47 100644
--- a/flow/_src/floww.i
+++ b/flow/_src/floww.i
@@ -51,7 +51,6 @@ threads="1"
%shared_ptr(flow::Initial);
%shared_ptr(flow::Dirichlet);
%shared_ptr(flow::Freestream);
-%shared_ptr(flow::Slip);
%shared_ptr(flow::Wake);
%shared_ptr(flow::Kutta);
%shared_ptr(flow::Blowing);
@@ -74,31 +73,6 @@ threads="1"
%include "wWakeElement.h"
%include "wKuttaElement.h"
-%feature("director") flow::F0PsC;
-%pythonappend flow::F0PsC::F0PsC "self.__disown__()" // only for directors
-%feature("director") flow::F0PsPhiInf;
-%pythonappend flow::F0PsPhiInf::F0PsPhiInf "self.__disown__()" // only for directors
-%nodefault flow::F0Ps; // is pure virtual but swig doesn't know it.
-%include "wF0Ps.h"
-
-%feature("director") flow::F0ElRho;
-%pythonappend flow::F0ElRho::F0ElRho "self.__disown__()" // only for directors
-%feature("director") flow::F0ElRhoL;
-%pythonappend flow::F0ElRhoL::F0ElRhoL "self.__disown__()" // only for directors
-%feature("director") flow::F0ElMach;
-%pythonappend flow::F0ElMach::F0ElMach "self.__disown__()" // only for directors
-%feature("director") flow::F0ElMachL;
-%pythonappend flow::F0ElMachL::F0ElMachL "self.__disown__()" // only for directors
-%feature("director") flow::F0ElCp;
-%pythonappend flow::F0ElCp::F0ElCp "self.__disown__()" // only for directors
-%feature("director") flow::F0ElCpL;
-%pythonappend flow::F0ElCpL::F0ElCpL "self.__disown__()" // only for directors
-%include "wF0El.h"
-
-%feature("director") flow::F1ElVi;
-%pythonappend flow::F1ElVi::F1ElVi "self.__disown__()" // only for directors
-%include "wF1El.h"
-
%immutable flow::Problem::msh; // read-only variable (no setter)
%include "wProblem.h"
diff --git a/flow/benchmark/onera.py b/flow/benchmark/onera.py
index 85aa46ab..9edc59a5 100644
--- a/flow/benchmark/onera.py
+++ b/flow/benchmark/onera.py
@@ -65,7 +65,6 @@ def main():
# define flow variables
alpha = 3.06*np.pi/180
M_inf = 0.839
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension
@@ -80,7 +79,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
diff --git a/flow/cases/coyote.py b/flow/cases/coyote.py
index cb311184..6994b90a 100644
--- a/flow/cases/coyote.py
+++ b/flow/cases/coyote.py
@@ -71,7 +71,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/n0012.py b/flow/cases/n0012.py
index a106ca95..52234831 100644
--- a/flow/cases/n0012.py
+++ b/flow/cases/n0012.py
@@ -50,7 +50,6 @@ def getParam():
p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
p['Max_it'] = 20 # Solver maximum number of iterations
p['Relaxation'] = 0.7 # Relaxation parameter for Picard solver
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/n0012_3.py b/flow/cases/n0012_3.py
index cc2c0344..c0d087e4 100644
--- a/flow/cases/n0012_3.py
+++ b/flow/cases/n0012_3.py
@@ -58,7 +58,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/n64A410.py b/flow/cases/n64A410.py
index e406b334..3b8fcaa7 100644
--- a/flow/cases/n64A410.py
+++ b/flow/cases/n64A410.py
@@ -52,7 +52,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/rae2822.py b/flow/cases/rae2822.py
index c5b9f0f2..a447e18d 100644
--- a/flow/cases/rae2822.py
+++ b/flow/cases/rae2822.py
@@ -55,7 +55,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/wbht.py b/flow/cases/wbht.py
index 69b2c484..2ab10ec0 100644
--- a/flow/cases/wbht.py
+++ b/flow/cases/wbht.py
@@ -83,7 +83,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/cases/wht.py b/flow/cases/wht.py
index 1efcd650..72deaf3e 100644
--- a/flow/cases/wht.py
+++ b/flow/cases/wht.py
@@ -74,7 +74,6 @@ def getParam():
p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
- p['M_crit'] = 5. # Critical Mach number above which density is damped
return p
def main():
diff --git a/flow/default.py b/flow/default.py
index a6b3e9b4..50b4bb56 100644
--- a/flow/default.py
+++ b/flow/default.py
@@ -42,29 +42,26 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
gmshWriter.save(msh.name)
return msh, gmshWriter
-def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False):
+def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False):
'''Initialize problem
'''
# create base formulation
pbl = flo.Problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref)
# add incompressible or compressible air medium
- if minf == 0:
- pbl.set(flo.Medium(msh, fld, flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha, beta)))
- else:
- pbl.set(flo.Medium(msh, fld, flo.F0ElRho(minf, mcrit), flo.F0ElMach(minf), flo.F0ElCp(minf), flo.F0PsPhiInf(dim, alpha, beta)))
+ pbl.set(flo.Medium(msh, fld, minf, dim, alpha, beta))
# add initial condition
- pbl.set(flo.Initial(msh, fld, flo.F0PsPhiInf(dim, alpha, beta)))
+ pbl.set(flo.Initial(msh, fld, dim, alpha, beta))
if dbc:
# Dirichlet BCs provided and explicitely imposed
- dirichlet = flo.Dirichlet(msh, upstr, flo.F0PsPhiInf(dim, alpha, beta))
+ dirichlet = flo.Dirichlet(msh, upstr, dim, alpha, beta)
pbl.add(dirichlet)
else:
# Dirichlet BCs not provided, Neumann imposed instead, and a DOF will be pinned automatically
dirichlet = None
- pbl.add(flo.Freestream(msh, upstr, flo.F1ElVi(dim, alpha, beta)))
+ pbl.add(flo.Freestream(msh, upstr, dim, alpha, beta))
# add freestream (Neumann) conditions
- pbl.add(flo.Freestream(msh, ff, flo.F1ElVi(dim, alpha, beta)))
- pbl.add(flo.Freestream(msh, dnstr, flo.F1ElVi(dim, alpha, beta)))
+ pbl.add(flo.Freestream(msh, ff, dim, alpha, beta))
+ pbl.add(flo.Freestream(msh, dnstr, dim, alpha, beta))
# add wake and Kutta conditions
if dim == 2 and te is not None:
wake = flo.Wake(msh, [wk, wk+'_', fld])
diff --git a/flow/scripts/config.py b/flow/scripts/config.py
index 8625c7a4..93c79c75 100644
--- a/flow/scripts/config.py
+++ b/flow/scripts/config.py
@@ -19,6 +19,7 @@
## Base class for configuration
# Adrien Crovato
+import math
import flow as f
import tbox
import tbox.gmsh as gmsh
@@ -39,11 +40,16 @@ class Config:
if p['Dim'] == 2:
if 'AoS' in p and p['AoS'] != 0:
raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
+ beta = 0
if 'Symmetry' in p:
raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
else:
- if 'AoS' in p and p['AoS'] != 0 and 'Symmetry' in p:
- raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
+ if 'AoS' in p:
+ if p['AoS'] != 0 and 'Symmetry' in p:
+ raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
+ beta = p['AoS']*math.pi/180
+ else:
+ beta = 0
# basic config
if p['Format'] == 'vtk':
try:
@@ -82,19 +88,14 @@ class Config:
print('\n')
# initialize the problem
- self.pbl = f.Problem(self.msh, p['Dim'], 0., 0., p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
- self.phiInfFun = f.F0PsPhiInf(p['Dim'], 0.)
- self.velInfFun = f.F1ElVi(p['Dim'], 0.)
+ self.pbl = f.Problem(self.msh, p['Dim'], 0., beta, p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
# add a medium "air"
- if p['M_inf'] == 0:
- self.pbl.set(f.Medium(self.msh, p['Fluid'], f.F0ElRhoL(), f.F0ElMachL(), f.F0ElCpL(), self.phiInfFun))
- else:
- self.pbl.set(f.Medium(self.msh, p['Fluid'], f.F0ElRho(p['M_inf'], p['M_crit']), f.F0ElMach(p['M_inf']), f.F0ElCp(p['M_inf']), self.phiInfFun))
+ self.pbl.set(f.Medium(self.msh, p['Fluid'], p['M_inf'], p['Dim'], 0., beta))
# add initial condition
- self.pbl.set(f.Initial(self.msh, p['Fluid'], self.phiInfFun))
+ self.pbl.set(f.Initial(self.msh, p['Fluid'], p['Dim'], 0., beta))
# add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
for i in range (0, len(p['Farfield'])):
- self.pbl.add(f.Freestream(self.msh, p['Farfield'][i], self.velInfFun))
+ self.pbl.add(f.Freestream(self.msh, p['Farfield'][i], p['Dim'], 0., beta))
# add solid boundaries
if p['Dim'] == 2:
self.bnd = f.Boundary(self.msh, [p['Wing'], p['Fluid']])
diff --git a/flow/scripts/polar.py b/flow/scripts/polar.py
index 7cb623ee..8d3e6942 100644
--- a/flow/scripts/polar.py
+++ b/flow/scripts/polar.py
@@ -40,10 +40,6 @@ class Polar(Config):
self.tag = None
self.mach = p['M_inf']
self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
- if 'AoS' in p:
- self.beta = p['AoS']*math.pi/180
- else:
- self.beta = 0
# basic init
Config.__init__(self, p)
@@ -58,10 +54,7 @@ class Polar(Config):
# define current angle of attack
alpha = alpha*math.pi/180
# update problem
- self.pbl.alpha = alpha
- self.pbl.beta = self.beta
- self.phiInfFun.update(alpha, self.beta)
- self.velInfFun.update(alpha, self.beta)
+ self.pbl.update(alpha)
# run the solver and save the results
print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
self.tms["solver"].start()
diff --git a/flow/scripts/trim.py b/flow/scripts/trim.py
index a14bea0c..1c3b0f5f 100644
--- a/flow/scripts/trim.py
+++ b/flow/scripts/trim.py
@@ -41,10 +41,6 @@ class Trim(Config):
self.mach = p['M_inf']
self.ClTgt = p['CL']
self.alpha = p['AoA']*math.pi/180
- if 'AoS' in p:
- self.beta = p['AoS']*math.pi/180
- else:
- self.beta = 0
self.dCl = p['dCL']
# basic init
Config.__init__(self, p)
@@ -56,10 +52,7 @@ class Trim(Config):
# define current angle of attack
print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
# update problem
- self.pbl.alpha = self.alpha
- self.pbl.beta = self.beta
- self.phiInfFun.update(self.alpha, self.beta)
- self.velInfFun.update(self.alpha, self.beta)
+ self.pbl.update(self.alpha)
# run the solver
self.tms["solver"].start()
success = self.solver.run()
diff --git a/flow/src/flow.h b/flow/src/flow.h
index 351ca549..27277a2b 100644
--- a/flow/src/flow.h
+++ b/flow/src/flow.h
@@ -36,7 +36,6 @@
*/
namespace flow
{
-
// Problem and B.C. handling
class Problem;
class Medium;
@@ -69,7 +68,19 @@ class Adjoint;
// General
class F0Ps;
class F0El;
+class F0ElC;
+class F0ElRho;
+class F0ElRhoL;
+class F0ElMach;
+class F0ElMachL;
+class F0ElCp;
+class F0ElCpL;
class F1El;
+class F1ElVi;
+class F1Ct;
+class F1CtDrag;
+class F1CtSide;
+class F1CtLift;
class Face;
class FaceEq;
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index 11eebf98..3eac871c 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -247,28 +247,22 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
// Multithread
tbb::spin_mutex mutex;
- // Reset gradients
+ // Reset and compute gradient of lift and drag
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) {
// Associated volume element
Element *eV = sur->svMap.at(e);
// Build
- Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, sol->pbl->medium->cP);
+ Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->medium->cP);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
- 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());
+ Eigen::RowVectorXd AeL = sol->pbl->dirL.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size());
+ Eigen::RowVectorXd AeD = sol->pbl->dirD.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size());
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
Node *nodj = eV->nodes[j];
@@ -307,26 +301,6 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
}
});
}
-
- // Wake contribution (through stabilisation term)
- // @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;
- 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)
- {
- Node *nodi = we->wkN[ii].second;
- dR(nodi->row) += bue(ii);
- dR(nodi->row) -= ble(ii);
- }
- });
- }
}
/**
@@ -335,17 +309,14 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
*/
void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
{
- // 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 (size_t i = 0; i < bnd->nodes.size(); ++i)
Cf += bnd->nLoads[i];
// Rotate to gradient of flow direction
- dD = Cf.dot(dirD);
- dL = Cf.dot(dirL);
+ dD = Cf.dot(sol->pbl->dirD.evalGrad());
+ dL = Cf.dot(sol->pbl->dirL.evalGrad());
}
/**
@@ -388,7 +359,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
}
}
// Assembly (supersonic)
- double mach = fluid->mach.eval(*e, sol->phi, 0);
+ double mach = fluid->mach->eval(*e, sol->phi, 0);
if (mach > sol->mCOv)
{
for (size_t i = 0; i < e->nodes.size(); ++i)
@@ -413,11 +384,10 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
// Wake
for (auto wake : sol->pbl->wBCs)
{
- Eigen::VectorXd vi = sol->pbl->computeVi();
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build elementary matrices
Eigen::MatrixXd Aeu, Ael, Aes;
- std::tie(Aeu, Ael, Aes) = WakeResidual::buildGradientMesh(*we, sol->phi, vi);
+ std::tie(Aeu, Ael, Aes) = WakeResidual::buildGradientMesh(*we, sol->phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < we->nRow; ++i)
@@ -511,15 +481,9 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd
// Multithread
tbb::spin_mutex mutex;
- // Reset gradients
+ // Reset and compute gradient of lift and drag
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 bnd : sol->pbl->bnds)
{
tbb::parallel_do(bnd->groups[0]->tag->elems.begin(), bnd->groups[0]->tag->elems.end(), [&](Element *e) {
@@ -527,16 +491,16 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd
Element *eV = bnd->svMap.at(e);
// Build
Eigen::MatrixXd Aes, Aev;
- std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, sol->pbl->medium->cP, sol->pbl->nDim);
+ 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 i = 0; i < e->nodes.size(); ++i)
{
// 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());
+ Eigen::RowVectorXd AesL = sol->pbl->dirL.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AesD = sol->pbl->dirD.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AevL = sol->pbl->dirL.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+ Eigen::RowVectorXd AevD = sol->pbl->dirD.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
for (size_t m = 0; m < sol->pbl->nDim; ++m)
{
// surface
diff --git a/flow/src/wAssign.cpp b/flow/src/wAssign.cpp
index 7056fcb6..82ecc405 100644
--- a/flow/src/wAssign.cpp
+++ b/flow/src/wAssign.cpp
@@ -23,30 +23,31 @@
using namespace tbox;
using namespace flow;
-Assign::Assign(std::shared_ptr<MshData> _msh, int no, F0Ps &_f) : Group(_msh, no), f(_f)
+Assign::Assign(std::shared_ptr<MshData> _msh, int no) : Group(_msh, no)
{
getNodes();
}
-Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name, F0Ps &_f) : Group(_msh, name), f(_f)
+Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name) : Group(_msh, name)
{
getNodes();
}
+Assign::~Assign()
+{
+ delete f;
+ std::cout << "~Assign()\n";
+}
/**
* @brief Prescribe value to all nodes (default implementation)
- *
- * Value can be constant or position dependent
- * @authors Adrien Crovato, Romain Boman
*/
void Assign::apply(std::vector<double> &vec)
{
for (auto n : nodes)
- vec[n->row] = f.eval(n->pos);
+ vec[n->row] = f->eval(n->pos);
}
/**
* @brief Get all nodes belonging to boundary tag
- * @authors Adrien Crovato, Romain Boman
*/
void Assign::getNodes()
{
@@ -58,31 +59,35 @@ void Assign::getNodes()
nodes.resize(std::distance(nodes.begin(), it));
}
-Initial::Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f) : Assign(_msh, no, _f)
+Initial::Initial(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta) : Assign(_msh, no)
{
+ f = new F0PsPhiInf(dim, alpha, beta);
}
-Initial::Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f)
+Initial::Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name)
{
+ f = new F0PsPhiInf(dim, alpha, beta);
}
void Initial::write(std::ostream &out) const
{
out << "flow::Initial condition on " << *tag << ")\n";
for (auto n : nodes)
- out << '\t' << *n << ": (val=" << f.eval(n->pos) << ")\n";
+ out << '\t' << *n << ": (val=" << f->eval(n->pos) << ")\n";
}
-Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f, bool pin) : Assign(_msh, no, _f)
+Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta, bool pin) : Assign(_msh, no)
{
+ f = new F0PsPhiInf(dim, alpha, beta);
// Only retain the first node, so that the DOF associated to this node will be pinned
if (pin)
this->nodes.resize(1);
}
-Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f)
+Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name)
{
+ f = new F0PsPhiInf(dim, alpha, beta);
}
void Dirichlet::write(std::ostream &out) const
{
out << "flow::Dirichlet boundary condition on " << *tag << ")\n";
for (auto n : nodes)
- out << '\t' << *n << ": (val=" << f.eval(n->pos) << ")\n";
+ out << '\t' << *n << ": (val=" << f->eval(n->pos) << ")\n";
}
\ No newline at end of file
diff --git a/flow/src/wAssign.h b/flow/src/wAssign.h
index e86d793b..784c51e9 100644
--- a/flow/src/wAssign.h
+++ b/flow/src/wAssign.h
@@ -36,11 +36,11 @@ class FLOW_API Assign : public tbox::Group
public:
std::vector<tbox::Node *> nodes; ///< nodes of the boundary
#ifndef SWIG
- F0Ps &f; ///< position-based function to compute the boundary value
+ F0Ps *f; ///< position-based function to compute the boundary value
- Assign(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f);
- Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
- virtual ~Assign() { std::cout << "~Assign()\n"; }
+ Assign(std::shared_ptr<tbox::MshData> _msh, int no);
+ Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name);
+ virtual ~Assign();
void apply(std::vector<double> &vec);
#endif
@@ -56,8 +56,8 @@ private:
class FLOW_API Initial : public Assign
{
public:
- Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f);
- Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
+ Initial(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0);
+ Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0);
virtual ~Initial() { std::cout << "~Initial()\n"; }
#ifndef SWIG
@@ -72,8 +72,8 @@ public:
class FLOW_API Dirichlet : public Assign
{
public:
- Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f, bool pin = false);
- Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
+ Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0, bool pin = false);
+ Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0);
virtual ~Dirichlet() { std::cout << "~Dirichlet()\n"; }
#ifndef SWIG
diff --git a/flow/src/wBlowing.cpp b/flow/src/wBlowing.cpp
index 29420b69..50992462 100644
--- a/flow/src/wBlowing.cpp
+++ b/flow/src/wBlowing.cpp
@@ -39,7 +39,6 @@ Blowing::~Blowing()
/**
* @brief Create data structure
- * @authors Adrien Crovato
*/
void Blowing::create()
{
@@ -55,14 +54,13 @@ void Blowing::create()
uE.resize(tag->elems.size());
for (size_t i = 0; i < tag->elems.size(); ++i)
{
- Fct0C *f = new Fct0C(0.);
- uE[i] = std::pair<Element *, Fct0C *>(tag->elems[i], f);
+ F0ElC *f = new F0ElC(0.);
+ uE[i] = std::pair<Element *, F0ElC *>(tag->elems[i], f);
}
}
/**
* @brief Set blowing velocity
- * @authors Adrien Crovato
*/
void Blowing::setU(size_t i, double ue)
{
diff --git a/flow/src/wBlowing.h b/flow/src/wBlowing.h
index de0f6400..c3706c31 100644
--- a/flow/src/wBlowing.h
+++ b/flow/src/wBlowing.h
@@ -19,7 +19,7 @@
#include "flow.h"
#include "wGroup.h"
-#include "wFct0.h"
+#include "wF0El.h"
#include <vector>
namespace flow
@@ -34,7 +34,7 @@ class FLOW_API Blowing : public tbox::Group
public:
std::vector<tbox::Node *> nodes; ///< nodes of Blowing
#ifndef SWIG
- std::vector<std::pair<tbox::Element *, tbox::Fct0C *>> uE; ///< blowing velocity for each element
+ std::vector<std::pair<tbox::Element *, F0ElC *>> uE; ///< blowing velocity for each element
#endif
Blowing(std::shared_ptr<tbox::MshData> _msh, int nos);
diff --git a/flow/src/wBlowingResidual.cpp b/flow/src/wBlowingResidual.cpp
index 15092468..7753c358 100644
--- a/flow/src/wBlowingResidual.cpp
+++ b/flow/src/wBlowingResidual.cpp
@@ -15,7 +15,7 @@
*/
#include "wBlowingResidual.h"
-#include "wFct0.h"
+#include "wF0El.h"
#include "wElement.h"
#include "wCache.h"
@@ -29,7 +29,7 @@ using namespace flow;
*
* b = \int psi * vn dV
*/
-Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, Fct0 const &f)
+Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f)
{
// Get pre-computed values
Cache &cache = e.getVCache();
diff --git a/flow/src/wBlowingResidual.h b/flow/src/wBlowingResidual.h
index 95e09a3e..8c77ef55 100644
--- a/flow/src/wBlowingResidual.h
+++ b/flow/src/wBlowingResidual.h
@@ -32,7 +32,7 @@ class FLOW_API BlowingResidual
{
public:
// Newton
- static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, tbox::Fct0 const &f);
+ static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, F0El const &f);
};
} // namespace flow
diff --git a/flow/src/wF0El.cpp b/flow/src/wF0El.cpp
index 56432195..0829796e 100644
--- a/flow/src/wF0El.cpp
+++ b/flow/src/wF0El.cpp
@@ -18,9 +18,30 @@
using namespace tbox;
using namespace flow;
+/**
+ * @brief Update the value
+ */
+void F0ElC::update(double _v)
+{
+ v = _v;
+}
+/**
+ * @brief Evaluate the constant
+ */
+double F0ElC::eval(Element const &e, std::vector<double> const &u, size_t k) const
+{
+ return v;
+}
+/**
+ * @brief Evaluate the gradient of the constant
+ */
+double F0ElC::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+{
+ return 0;
+}
+
/**
* @brief Evaluate the nonlinear density (constant over element)
- * @authors Adrien Crovato
*/
double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -41,9 +62,8 @@ double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) c
}
/**
* @brief Evaluate the nonlinear density derivative (constant over element)
- * @authors Adrien Crovato
*/
-double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRho::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
double gradU = e.computeGradient(u, k).norm(); // norm of velocity
@@ -60,14 +80,9 @@ double F0ElRho::evalD(Element const &e, std::vector<double> const &u, size_t k)
return -sqrt(a * a * a * a * a) / ((1 + beta * (gradU / gradUC - 1)) * (1 + beta * (gradU / gradUC - 1))) * beta / gradUC / gradU; // has to be multiplied by grad u to recover true derivative
}
}
-void F0ElRho::write(std::ostream &out) const
-{
- out << "F0ElRho";
-}
/**
* @brief Evaluate the linear density (constant over element)
- * @authors Adrien Crovato
*/
double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -75,20 +90,14 @@ double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k)
}
/**
* @brief Evaluate the linear density derivative (constant over element)
- * @authors Adrien Crovato
*/
-double F0ElRhoL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
return 0.;
}
-void F0ElRhoL::write(std::ostream &out) const
-{
- out << "F0ElRhoL";
-}
/**
* @brief Evaluate the nonlinear Mach number (constant over element)
- * @authors Adrien Crovato
*/
double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -100,9 +109,8 @@ double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k)
}
/**
* @brief Evaluate the nonlinear Mach number derivative (constant over element)
- * @authors Adrien Crovato
*/
-double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMach::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
double gradU2 = gradU.squaredNorm(); // velocity squared
@@ -110,14 +118,9 @@ double F0ElMach::evalD(Element const &e, std::vector<double> const &u, size_t k)
return 1 / sqrt(gradU2 * a2) + 0.5 * (gamma - 1) * sqrt(gradU2) / sqrt(a2 * a2 * a2); // has to be multiplied by grad u to recover true derivative
}
-void F0ElMach::write(std::ostream &out) const
-{
- out << "F0ElMach";
-}
/**
* @brief Evaluate the linear Mach number (constant over element)
- * @authors Adrien Crovato
*/
double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -125,20 +128,14 @@ double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k)
}
/**
* @brief Evaluate the linear Mach number derivative (constant over element)
- * @authors Adrien Crovato
*/
-double F0ElMachL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
return 0.;
}
-void F0ElMachL::write(std::ostream &out) const
-{
- out << "F0ElMachL";
-}
/**
* @brief Evaluate the nonlinear pressure coefficient (constant over element)
- * @authors Adrien Crovato
*/
double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -154,9 +151,8 @@ double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) co
}
/**
* @brief Evaluate the nonlinear pressure coefficient derivative (constant over element)
- * @authors Adrien Crovato
*/
-double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
double gradU2 = gradU.squaredNorm(); // velocity squared
@@ -168,14 +164,9 @@ double F0ElCp::evalD(Element const &e, std::vector<double> const &u, size_t k) c
else
return -2 * sqrt(a * a * a * a * a); // has to be multiplied by grad u to recover true derivative
}
-void F0ElCp::write(std::ostream &out) const
-{
- out << "F0ElCp";
-}
/**
* @brief Evaluate the linear pressure coefficient (constant over element)
- * @authors Adrien Crovato
*/
double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -183,11 +174,7 @@ double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) c
double gradU2 = gradU.squaredNorm(); // velocity squared
return 1 - gradU2;
}
-double F0ElCpL::evalD(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCpL::evalGrad(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
-}
-void F0ElCpL::write(std::ostream &out) const
-{
- out << "F0ElCpL";
}
\ No newline at end of file
diff --git a/flow/src/wF0El.h b/flow/src/wF0El.h
index a92e4d4b..d046d991 100644
--- a/flow/src/wF0El.h
+++ b/flow/src/wF0El.h
@@ -18,153 +18,149 @@
#define WF0EL_H
#include "flow.h"
-#include "wFct0.h"
+#include "wObserver.h"
#include "wElement.h"
#include <vector>
namespace flow
{
+/**
+ * @brief Scalar function to be integrated over an element
+ */
+class FLOW_API F0El : public fwk::Observer
+{
+public:
+ F0El() : Observer() {}
+ virtual ~F0El() {}
+
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const = 0;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const = 0;
+};
+
+/**
+ * @brief Constant function
+ */
+class FLOW_API F0ElC : public F0El
+{
+ double v;
+
+public:
+ F0ElC(double _v) : F0El(), v(_v) {}
+ virtual void update(double _v) override;
+
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
+};
+
/**
* @brief Nonlinear density
- *
+ *
* The density is limited (Padé approximation) for large value of velocity and
* particularized for diatomic gas. The input _mC2 is the square of the critical
* Mach number used to compute the critical velocity.
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElRho : public tbox::Fct0
+class FLOW_API F0ElRho : public F0El
{
double gamma; ///< heat capacity ratio (diatomic gas only)
double mInf; ///< freestream Mach number
double gradUC; ///< critical velocity
public:
- F0ElRho(double _mInf, double _mC2) : Fct0(), gamma(1.4), mInf(_mInf)
+ F0ElRho(double _mInf, double _mC2 = 5) : F0El(), gamma(1.4), mInf(_mInf)
{
gradUC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2));
}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
/**
* @brief Linear density (constant)
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElRhoL : public tbox::Fct0
+class FLOW_API F0ElRhoL : public F0El
{
public:
- F0ElRhoL() : Fct0() {}
+ F0ElRhoL() : F0El() {}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
/**
* @brief Nonlinear Mach number
- *
+ *
* Particularized for diatomic gas
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElMach : public tbox::Fct0
+class FLOW_API F0ElMach : public F0El
{
double gamma; ///< heat capacity ratio (diatomic gas only)
double mInf; ///< freestream Mach number
public:
- F0ElMach(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {}
+ F0ElMach(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
/**
* @brief Linear Mach number (constant)
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElMachL : public tbox::Fct0
+class FLOW_API F0ElMachL : public F0El
{
public:
- F0ElMachL() : Fct0() {}
+ F0ElMachL() : F0El() {}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
/**
* @brief Nonlinear pressure coefficient
- *
+ *
* Particularized for diatomic gas
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElCp : public tbox::Fct0
+class FLOW_API F0ElCp : public F0El
{
double gamma; ///< heat capacity ratio (diatomic gas only)
double mInf; ///< freestream Mach number
public:
- F0ElCp(double _mInf) : Fct0(), gamma(1.4), mInf(_mInf) {}
+ F0ElCp(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
/**
* @brief Linear pressure coefficient
- * @authors Adrien Crovato
*/
-class FLOW_API F0ElCpL : public tbox::Fct0
+class FLOW_API F0ElCpL : public F0El
{
public:
- F0ElCpL() : Fct0() {}
+ F0ElCpL() : F0El() {}
-#ifndef SWIG
- virtual double eval(tbox::Element const &e,
- std::vector<double> const &u,
+ virtual double eval(tbox::Element const &e, std::vector<double> const &u,
size_t k) const override;
- 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;
-#endif
+ virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
} // namespace flow
diff --git a/flow/src/wF0Ps.cpp b/flow/src/wF0Ps.cpp
index bc25ca8c..7c8a0425 100644
--- a/flow/src/wF0Ps.cpp
+++ b/flow/src/wF0Ps.cpp
@@ -15,23 +15,9 @@
*/
#include "wF0Ps.h"
-#include <sstream>
using namespace tbox;
using namespace flow;
-/**
- * @brief Return constant function
- * @authors Adrien Crovato
- */
-double F0PsC::eval(Eigen::Vector3d const &pos) const
-{
- return v;
-}
-void F0PsC::write(std::ostream &out) const
-{
- out << "F0PsC";
-}
-
F0PsPhiInf::F0PsPhiInf(int _nDim, double _alpha, double _beta) : F0Ps(),
nDim(_nDim),
alpha(_alpha),
@@ -48,17 +34,14 @@ F0PsPhiInf::F0PsPhiInf(int _nDim, double _alpha, double _beta) : F0Ps(),
throw std::runtime_error("F0PsPhiInf::beta can be nonzero only for 3D cases!\n");
}
/**
- * @brief Update the angle of attack and the angle of sideslip
- * @authors Adrien Crovato
+ * @brief Update the angle of attack
*/
-void F0PsPhiInf::update(double _alpha, double _beta)
+void F0PsPhiInf::update(double _alpha)
{
alpha = _alpha;
- beta = _beta;
}
/**
* @brief Evaluate the freestream potential
- * @authors Adrien Crovato
*/
double F0PsPhiInf::eval(Eigen::Vector3d const &pos) const
{
@@ -67,7 +50,3 @@ double F0PsPhiInf::eval(Eigen::Vector3d const &pos) const
else
return cos(alpha) * cos(beta) * pos(0) + sin(beta) * pos(1) + sin(alpha) * cos(beta) * pos(2);
}
-void F0PsPhiInf::write(std::ostream &out) const
-{
- out << "F0PsPhiInf";
-}
diff --git a/flow/src/wF0Ps.h b/flow/src/wF0Ps.h
index b8b17607..b68f2dd8 100644
--- a/flow/src/wF0Ps.h
+++ b/flow/src/wF0Ps.h
@@ -18,7 +18,7 @@
#define WF0PS_H
#include "flow.h"
-#include "wObject.h"
+#include "wObserver.h"
#include <Eigen/Dense>
namespace flow
@@ -26,44 +26,18 @@ namespace flow
/**
* @brief Scalar functions depending on spatial position
- * @authors Adrien Crovato
*/
-class FLOW_API F0Ps : public fwk::wObject
+class FLOW_API F0Ps : public fwk::Observer
{
public:
-#ifndef SWIG
- F0Ps() : wObject()
- {
- }
-#endif
- virtual ~F0Ps()
- {
- }
-#ifndef SWIG
- virtual double eval(Eigen::Vector3d const &pos) const = 0;
-#endif
-};
-
-/**
- * @brief Constant function
- * @authors Adrien Crovato
- */
-class FLOW_API F0PsC : public F0Ps
-{
- double v;
-
-public:
- F0PsC(double _v) : F0Ps(), v(_v) {}
+ F0Ps() : Observer() {}
+ virtual ~F0Ps() {}
-#ifndef SWIG
- virtual double eval(Eigen::Vector3d const &pos) const override;
- virtual void write(std::ostream &out) const override;
-#endif
+ virtual double eval(Eigen::Vector3d const &pos) const = 0;
};
/**
* @brief Freestream potential
- * @authors Adrien Crovato
*/
class FLOW_API F0PsPhiInf : public F0Ps
{
@@ -73,12 +47,9 @@ class FLOW_API F0PsPhiInf : public F0Ps
public:
F0PsPhiInf(int _nDim, double _alpha, double _beta = 0);
- void update(double _alpha, double _beta = 0);
+ virtual void update(double _alpha) override;
-#ifndef SWIG
virtual double eval(Eigen::Vector3d const &pos) const override;
- virtual void write(std::ostream &out) const override;
-#endif
};
} // namespace flow
diff --git a/flow/src/wF1Ct.cpp b/flow/src/wF1Ct.cpp
new file mode 100644
index 00000000..57769dfe
--- /dev/null
+++ b/flow/src/wF1Ct.cpp
@@ -0,0 +1,154 @@
+/*
+ * 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 "wF1Ct.h"
+using namespace tbox;
+using namespace flow;
+
+F1CtDrag::F1CtDrag(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta)
+{
+ // Sanity check
+ if (nDim != 2 && nDim != 3)
+ {
+ std::stringstream err;
+ err << "F1CtDrag::nDim should be equal to 2 or 3 but " << nDim << " was given!\n";
+ throw std::runtime_error(err.str());
+ }
+ if (nDim != 3 && beta != 0)
+ throw std::runtime_error("F1CtDrag::beta can be nonzero only for 3D cases!\n");
+ // Set the velocity
+ this->update(alpha);
+}
+/**
+ * @brief Update the direction and its gradient
+ */
+void F1CtDrag::update(double alpha)
+{
+ if (nDim == 2)
+ {
+ v = Eigen::Vector3d(cos(alpha), sin(alpha), 0) / sRef;
+ dv = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / sRef;
+ }
+ else
+ {
+ v = Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)) / sRef;
+ dv = Eigen::Vector3d(-sin(alpha) * cos(beta), sin(beta), cos(alpha) * cos(beta)) / sRef;
+ }
+}
+/**
+ * @brief Return the direction
+ */
+Eigen::Vector3d F1CtDrag::eval() const
+{
+ return v;
+}
+/**
+ * @brief Return the gradient of the direction with respect to angle of attack
+ */
+Eigen::Vector3d F1CtDrag::evalGrad() const
+{
+ return dv;
+}
+
+F1CtSide::F1CtSide(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta)
+{
+ // Sanity check
+ if (nDim != 2 && nDim != 3)
+ {
+ std::stringstream err;
+ err << "F1CtSide::nDim should be equal to 2 or 3 but " << nDim << " was given!\n";
+ throw std::runtime_error(err.str());
+ }
+ if (nDim != 3 && beta != 0)
+ throw std::runtime_error("F1CtSide::beta can be nonzero only for 3D cases!\n");
+ // Set the velocity
+ this->update(alpha);
+}
+/**
+ * @brief Update the direction and its gradient
+ */
+void F1CtSide::update(double alpha)
+{
+ if (nDim == 2)
+ {
+ v = Eigen::Vector3d(0, 0, 0);
+ dv = Eigen::Vector3d(0, 0, 0);
+ }
+ else
+ {
+ v = Eigen::Vector3d(-cos(alpha) * sin(beta), cos(beta), -sin(alpha) * sin(beta)) / sRef;
+ dv = Eigen::Vector3d(sin(alpha) * sin(beta), cos(beta), -cos(alpha) * sin(beta)) / sRef;
+ }
+}
+/**
+ * @brief Return the direction
+ */
+Eigen::Vector3d F1CtSide::eval() const
+{
+ return v;
+}
+/**
+ * @brief Return the gradient of the direction with respect to angle of attack
+ */
+Eigen::Vector3d F1CtSide::evalGrad() const
+{
+ return dv;
+}
+
+F1CtLift::F1CtLift(int _nDim, double _sRef, double alpha, double _beta) : F1Ct(), nDim(_nDim), sRef(_sRef), beta(_beta)
+{
+ // Sanity check
+ if (nDim != 2 && nDim != 3)
+ {
+ std::stringstream err;
+ err << "F1CtLift::nDim should be equal to 2 or 3 but " << nDim << " was given!\n";
+ throw std::runtime_error(err.str());
+ }
+ if (nDim != 3 && beta != 0)
+ throw std::runtime_error("F1CtLift::beta can be nonzero only for 3D cases!\n");
+ // Set the velocity
+ this->update(alpha);
+}
+/**
+ * @brief Update the direction and its gradient
+ */
+void F1CtLift::update(double alpha)
+{
+ if (nDim == 2)
+ {
+ v = Eigen::Vector3d(-sin(alpha), cos(alpha), 0) / sRef;
+ dv = Eigen::Vector3d(-cos(alpha), -sin(alpha), 0) / sRef;
+ }
+ else
+ {
+ v = Eigen::Vector3d(-sin(alpha), 0, cos(alpha)) / sRef;
+ dv = Eigen::Vector3d(-cos(alpha), 0, -sin(alpha)) / sRef;
+ }
+}
+/**
+ * @brief Return the direction
+ */
+Eigen::Vector3d F1CtLift::eval() const
+{
+ return v;
+}
+/**
+ * @brief Return the gradient of the direction with respect to angle of attack
+ */
+Eigen::Vector3d F1CtLift::evalGrad() const
+{
+ return dv;
+}
diff --git a/flow/src/wF1Ct.h b/flow/src/wF1Ct.h
new file mode 100644
index 00000000..7af9fe21
--- /dev/null
+++ b/flow/src/wF1Ct.h
@@ -0,0 +1,99 @@
+/*
+ * 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 WF1CT_H
+#define WF1CT_H
+
+#include "flow.h"
+#include "wObserver.h"
+#include <Eigen/Dense>
+
+namespace flow
+{
+
+/**
+ * @brief Constant vector function
+ */
+class FLOW_API F1Ct : public fwk::Observer
+{
+public:
+ F1Ct() : Observer() {}
+ virtual ~F1Ct() {}
+
+ virtual Eigen::Vector3d eval() const = 0;
+ virtual Eigen::Vector3d evalGrad() const = 0;
+};
+
+/**
+ * @brief Drag direction
+ */
+class FLOW_API F1CtDrag : public F1Ct
+{
+ int nDim; ///< number of dimensions
+ double sRef; ///< reference area
+ double beta; ///< angle of sideslip
+ Eigen::Vector3d v; ///< drag direction unit vector
+ Eigen::Vector3d dv; ///< gradient of drag direction unit vector (wrt alpha)
+
+public:
+ F1CtDrag(int _nDim, double _sRef, double alpha, double _beta = 0);
+ void update(double alpha);
+
+ virtual Eigen::Vector3d eval() const override;
+ virtual Eigen::Vector3d evalGrad() const override;
+};
+
+/**
+ * @brief Sideforce direction
+ */
+class FLOW_API F1CtSide : public F1Ct
+{
+ int nDim; ///< number of dimensions
+ double sRef; ///< reference area
+ double beta; ///< angle of sideslip
+ Eigen::Vector3d v; ///< sideforce direction unit vector
+ Eigen::Vector3d dv; ///< gradient of sideforce direction unit vector (wrt alpha)
+
+public:
+ F1CtSide(int _nDim, double _sRef, double alpha, double _beta = 0);
+ void update(double alpha);
+
+ virtual Eigen::Vector3d eval() const override;
+ virtual Eigen::Vector3d evalGrad() const override;
+};
+
+/**
+ * @brief Lift direction
+ */
+class FLOW_API F1CtLift : public F1Ct
+{
+ int nDim; ///< number of dimensions
+ double sRef; ///< reference area
+ double beta; ///< angle of sideslip
+ Eigen::Vector3d v; ///< lift direction unit vector
+ Eigen::Vector3d dv; ///< gradient of lift direction unit vector (wrt alpha)
+
+public:
+ F1CtLift(int _nDim, double _sRef, double alpha, double _beta = 0);
+ void update(double alpha);
+
+ virtual Eigen::Vector3d eval() const override;
+ virtual Eigen::Vector3d evalGrad() const override;
+};
+
+} // namespace flow
+
+#endif //WF1CT_H
diff --git a/flow/src/wF1El.cpp b/flow/src/wF1El.cpp
index a65f35c8..d3d41fb7 100644
--- a/flow/src/wF1El.cpp
+++ b/flow/src/wF1El.cpp
@@ -18,7 +18,7 @@
using namespace tbox;
using namespace flow;
-F1ElVi::F1ElVi(int _nDim, double alpha, double beta) : Fct1(), nDim(_nDim)
+F1ElVi::F1ElVi(int _nDim, double alpha, double _beta) : F1El(), nDim(_nDim), beta(_beta)
{
// Sanity check
if (nDim != 2 && nDim != 3)
@@ -30,13 +30,12 @@ F1ElVi::F1ElVi(int _nDim, double alpha, double beta) : Fct1(), nDim(_nDim)
if (nDim != 3 && beta != 0)
throw std::runtime_error("F1ElVi::beta can be nonzero only for 3D cases!\n");
// Set the velocity
- this->update(alpha, beta);
+ this->update(alpha);
}
/**
* @brief Update the velocity vector
- * @authors Adrien Crovato
*/
-void F1ElVi::update(double alpha, double beta)
+void F1ElVi::update(double alpha)
{
if (nDim == 2)
{
@@ -51,7 +50,6 @@ void F1ElVi::update(double alpha, double beta)
}
/**
* @brief Return the freestream velocity
- * @authors Adrien Crovato
*/
Eigen::Vector3d F1ElVi::eval(Element const &e, std::vector<double> const &u, size_t k) const
{
@@ -60,11 +58,7 @@ Eigen::Vector3d F1ElVi::eval(Element const &e, std::vector<double> const &u, siz
/**
* @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
+Eigen::Vector3d F1ElVi::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
{
return dv;
}
-void F1ElVi::write(std::ostream &out) const
-{
- out << "F1ElVi";
-}
\ No newline at end of file
diff --git a/flow/src/wF1El.h b/flow/src/wF1El.h
index 7d869712..482b13e3 100644
--- a/flow/src/wF1El.h
+++ b/flow/src/wF1El.h
@@ -18,36 +18,46 @@
#define WF1EL_H
#include "flow.h"
-#include "wFct1.h"
+#include "wObserver.h"
#include <vector>
#include <Eigen/Dense>
namespace flow
{
+/**
+ * @brief Vector function to be integrated over an element
+ */
+class FLOW_API F1El : public fwk::Observer
+{
+public:
+ F1El() : Observer() {}
+ virtual ~F1El() {}
+
+ virtual Eigen::Vector3d eval(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const = 0;
+ virtual Eigen::Vector3d evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const = 0;
+};
+
/**
* @brief Freestream velocity
- * @authors Adrien Crovato
*/
-class FLOW_API F1ElVi : public tbox::Fct1
+class FLOW_API F1ElVi : public F1El
{
int nDim; ///< number of dimensions
+ double beta; ///< angle of sideslip
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);
+ F1ElVi(int _nDim, double alpha, double _beta = 0);
+ virtual void update(double alpha) override;
-#ifndef SWIG
- virtual Eigen::Vector3d eval(tbox::Element const &e,
- std::vector<double> const &u,
+ 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
+ virtual Eigen::Vector3d evalGrad(tbox::Element const &e, std::vector<double> const &u,
+ size_t k) const override;
};
} // namespace flow
diff --git a/flow/src/wFreestream.cpp b/flow/src/wFreestream.cpp
index b319b5fe..3ade264b 100644
--- a/flow/src/wFreestream.cpp
+++ b/flow/src/wFreestream.cpp
@@ -15,17 +15,23 @@
*/
#include "wFreestream.h"
-#include "wProblem.h"
#include "wTag.h"
using namespace tbox;
using namespace flow;
-Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, Fct1 &_f) : Group(_msh, no), f(_f)
+Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, int dim, double alpha, double beta) : Group(_msh, no)
{
+ f = new F1ElVi(dim, alpha, beta);
}
-Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, Fct1 &_f) : Group(_msh, name), f(_f)
+Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Group(_msh, name)
{
+ f = new F1ElVi(dim, alpha, beta);
+}
+Freestream::~Freestream()
+{
+ delete f;
+ std::cout << "~Freestream()\n";
}
void Freestream::write(std::ostream &out) const
diff --git a/flow/src/wFreestream.h b/flow/src/wFreestream.h
index 1725a254..96a78719 100644
--- a/flow/src/wFreestream.h
+++ b/flow/src/wFreestream.h
@@ -19,24 +19,24 @@
#include "flow.h"
#include "wGroup.h"
-#include "wFct1.h"
+#include "wF1EL.h"
namespace flow
{
/**
- * @brief Manage freestream (Freestream) boundary condtion
+ * @brief Manage freestream (Neumann) boundary condtion
* @authors Adrien Crovato
*/
class FLOW_API Freestream : public tbox::Group
{
public:
#ifndef SWIG
- tbox::Fct1 &f; ///< vector valued function to compute flux
+ F1ElVi *f; ///< vector valued function to compute flux
#endif
- Freestream(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct1 &_f);
- Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct1 &_f);
- virtual ~Freestream() { std::cout << "~Freestream()\n"; }
+ Freestream(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta = 0.0);
+ Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta = 0.0);
+ virtual ~Freestream();
#ifndef SWIG
virtual void write(std::ostream &out) const override;
diff --git a/flow/src/wFreestreamResidual.cpp b/flow/src/wFreestreamResidual.cpp
index b2678b45..04c8a247 100644
--- a/flow/src/wFreestreamResidual.cpp
+++ b/flow/src/wFreestreamResidual.cpp
@@ -36,7 +36,7 @@ Eigen::VectorXd FreestreamResidual::build(Element const &e, std::vector<double>
Gauss &gauss = cache.getVGauss();
// Normal velocity
- double vn = fs.f.eval(e, phi, 0).dot(e.getNormal());
+ double vn = fs.f->eval(e, phi, 0).dot(e.getNormal());
// Integral of mass flux (density times normal velocity) through freestream boundary
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
@@ -57,7 +57,7 @@ Eigen::VectorXd FreestreamResidual::buildGradientAoa(tbox::Element const &e, std
Gauss &gauss = cache.getVGauss();
// Gradient of normal velocity
- double dvn = fs.f.evalD(e, phi, 0).dot(e.getNormal());
+ double dvn = fs.f->evalGrad(e, phi, 0).dot(e.getNormal());
// Build mass flux gradient
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index 1113e3d8..401de65a 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -15,7 +15,7 @@
*/
#include "wLoadFunctional.h"
-#include "wFct0.h"
+#include "wF0El.h"
#include "wElement.h"
#include "wCache.h"
@@ -29,7 +29,7 @@ using namespace flow;
*
* b = \int psi * cp * n dS
*/
-Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp)
+Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp)
{
// Get pre-computed values
Cache &cache = e.getVCache();
@@ -55,14 +55,14 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::
*
* A = \int psi * dcp * n dS
*/
-Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element const &eV, std::vector<double> const &phi, Fct0 const &cp)
+Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp)
{
// Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
// Gradient of pressure force coefficient in associated volume
- Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid
+ Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalGrad(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid
// Integrate cp on the surface
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
@@ -84,7 +84,7 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
* + \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)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(Element const &e, Element const &eV, std::vector<double> const &phi, F0El const &cp, int nDim)
{
// Get pre-computed values
Cache &cache = e.getVCache();
@@ -94,7 +94,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
double cP = -cp.eval(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid)
Eigen::Vector3d nrm = e.getNormal(); // unit normal
// Gradients
- double dCp = -cp.evalD(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid)
+ double dCp = -cp.evalGrad(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid)
Eigen::VectorXd dPhi = eV.computeGradient(phi, 0); // potential
Eigen::MatrixXd const &iJ = eV.getJinv(0); // inverse Jacobian
std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian
diff --git a/flow/src/wLoadFunctional.h b/flow/src/wLoadFunctional.h
index 4a1a25a8..105fe2b7 100644
--- a/flow/src/wLoadFunctional.h
+++ b/flow/src/wLoadFunctional.h
@@ -32,10 +32,10 @@ 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);
+ static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El 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);
+ static Eigen::MatrixXd buildGradientFlow(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El const &cp);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eV, std::vector<double> const &phi, F0El const &cp, int nDim);
};
} // namespace flow
diff --git a/flow/src/wMedium.cpp b/flow/src/wMedium.cpp
index 9d4fecc2..c03d9a81 100644
--- a/flow/src/wMedium.cpp
+++ b/flow/src/wMedium.cpp
@@ -21,27 +21,54 @@
using namespace tbox;
using namespace flow;
-Medium::Medium(std::shared_ptr<MshData> _msh, int no, Fct0 &_rho, Fct0 &_mach,
- Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, no), rho(_rho),
- mach(_mach), cP(_cP), phiInf(_phiInf)
+Medium::Medium(std::shared_ptr<MshData> _msh, int no,
+ double mInf, int dim, double alpha, double beta) : Group(_msh, no)
{
+ initFun(mInf, dim, alpha, beta);
createMap();
std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes."
<< std::endl
<< std::endl;
}
Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name,
- Fct0 &_rho, Fct0 &_mach, Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, name), rho(_rho), mach(_mach), cP(_cP), phiInf(_phiInf)
+ double mInf, int dim, double alpha, double beta) : Group(_msh, name)
{
+ initFun(mInf, dim, alpha, beta);
createMap();
std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes."
<< std::endl
<< std::endl;
}
+Medium::~Medium()
+{
+ delete rho;
+ delete mach;
+ delete cP;
+ std::cout << "~Medium()\n";
+}
+
+/**
+ * Initialize functions
+ */
+void Medium::initFun(double mInf, int dim, double alpha, double beta)
+{
+ if (mInf == 0)
+ {
+ rho = new F0ElRhoL();
+ mach = new F0ElMachL();
+ cP = new F0ElCpL();
+ }
+ else
+ {
+ rho = new F0ElRho(mInf);
+ mach = new F0ElMach(mInf);
+ cP = new F0ElCp(mInf);
+ }
+ phiInf = new F0PsPhiInf(dim, alpha, beta);
+}
/**
* @brief Create data structure
- * @authors Adrien Crovato
*/
void Medium::createMap()
{
diff --git a/flow/src/wMedium.h b/flow/src/wMedium.h
index f84b815d..7bfbe691 100644
--- a/flow/src/wMedium.h
+++ b/flow/src/wMedium.h
@@ -19,8 +19,9 @@
#include "flow.h"
#include "wGroup.h"
-#include "wFct0.h"
+#include "wF0El.h"
#include "wF0Ps.h"
+#include <map>
namespace flow
{
@@ -31,24 +32,26 @@ namespace flow
*/
class FLOW_API Medium : public tbox::Group
{
+private:
+ size_t nCnt; ///< number of nodes
public:
std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< map each elements connected to a node
- size_t nCnt; ///< number of nodes
#ifndef SWIG
- tbox::Fct0 ρ ///< density
- tbox::Fct0 &mach; ///< Mach number
- tbox::Fct0 &cP; ///< pressure coefficient
- flow::F0PsPhiInf &phiInf; ///< freestream potential
+ F0El *rho; ///< density
+ F0El *mach; ///< Mach number
+ F0El *cP; ///< pressure coefficient
+ F0PsPhiInf *phiInf; ///< freestream potential
#endif
std::vector<std::pair<tbox::Element *, std::vector<tbox::Element *>>> adjMap; ///< map of adjacent elements
- Medium(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf);
- Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf);
- virtual ~Medium() { std::cout << "~Medium()\n"; }
+ Medium(std::shared_ptr<tbox::MshData> _msh, int no, double mInf, int dim, double alpha, double beta);
+ Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, double mInf, int dim, double alpha, double beta);
+ virtual ~Medium();
#ifndef SWIG
+ void initFun(double mInf, int dim, double alpha, double beta);
void createMap();
virtual void write(std::ostream &out) const override;
#endif
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index ec0a579a..4bff9008 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -287,7 +287,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
}
}
// Assembly (supersonic)
- double mach = fluid->mach.eval(*e, phi, 0);
+ double mach = fluid->mach->eval(*e, phi, 0);
if (mach > mCOv)
{
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
@@ -309,11 +309,10 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
tms["1-Jkutta"].start();
for (auto wake : pbl->wBCs)
{
- Eigen::VectorXd vi = pbl->computeVi();
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build Ku, Kl
Eigen::MatrixXd Aue, Ale;
- std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi, vi);
+ std::tie(Aue, Ale) = WakeResidual::buildGradientFlow(*we, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -431,7 +430,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
// Apply blowing (transpiration) BCs
for (auto bBC : pbl->bBCs)
{
- tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) {
+ tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
// Build elementary flux vector
Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second);
@@ -451,11 +450,10 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
// Nishida Thesis, pp. 29-30, Eq. 2.18
for (auto wake : pbl->wBCs)
{
- Eigen::VectorXd vi = pbl->computeVi();
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build bu, bl
Eigen::VectorXd bue, ble;
- std::tie(bue, ble) = WakeResidual::build(*we, phi, vi);
+ std::tie(bue, ble) = WakeResidual::build(*we, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
diff --git a/flow/src/wPicard.cpp b/flow/src/wPicard.cpp
index 8b66292b..4a3a890b 100644
--- a/flow/src/wPicard.cpp
+++ b/flow/src/wPicard.cpp
@@ -226,11 +226,10 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
tms["1-Akutta"].start();
for (auto wake : pbl->wBCs)
{
- Eigen::VectorXd vi = pbl->computeVi();
tbb::parallel_do(wake->wEle.begin(), wake->wEle.end(), [&](WakeElement *we) {
// Build Ku, Kl
Eigen::MatrixXd Aue, Ale;
- std::tie(Aue, Ale) = WakeResidual::buildFixed(*we, phi, vi);
+ std::tie(Aue, Ale) = WakeResidual::buildFixed(*we, phi);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < we->nRow; ++ii)
@@ -293,7 +292,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
// Apply blowing (transpiration) BCs
for (auto bBC : pbl->bBCs)
{
- tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, Fct0C *> p) {
+ tbb::parallel_do(bBC->uE.begin(), bBC->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
// Build elementary flux vector
Eigen::VectorXd be = BlowingResidual::build(*p.first, phi, *p.second);
diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp
index 1e683747..9122b663 100644
--- a/flow/src/wPotentialResidual.cpp
+++ b/flow/src/wPotentialResidual.cpp
@@ -42,7 +42,7 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
// Shape functions gradient
Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
// Elementary stiffness matrix
- A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k);
+ A += dSf.transpose() * dSf * fluid.rho->eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k);
}
return A;
}
@@ -51,7 +51,6 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
* @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)
{
@@ -62,14 +61,14 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
// 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) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
+ b += fluid.rho->eval(e, phi, k) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
// Supersonic contribution
- double mach = fluid.mach.eval(e, phi, 0);
+ 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
+ double rhoU = fluid.rho->eval(eU, phi, 0); // density on upwind element
// scale 1st term
b *= 1 - mu;
// contribution of upwind element
@@ -87,7 +86,6 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
* + \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)
{
@@ -106,13 +104,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
Eigen::VectorXd dPhi = e.computeGradient(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;
+ A1 += fluid.rho->eval(e, phi, k) * dSf.transpose() * dSf * wdj;
+ A1 += fluid.rho->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
}
// Supersonic contribution
Eigen::MatrixXd A2;
- double mach = fluid.mach.eval(e, phi, 0);
+ double mach = fluid.mach->eval(e, phi, 0);
if (mach > mCO)
{
// switching function and gradient
@@ -120,8 +118,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
// density and gradients on upwind elements
Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
- double rhoU = fluid.rho.eval(eU, phi, 0);
- double dRhoU = fluid.rho.evalD(eU, phi, 0);
+ double rhoU = fluid.rho->eval(eU, phi, 0);
+ double dRhoU = fluid.rho->evalGrad(eU, phi, 0);
Eigen::MatrixXd const &dSfU = eU.getJinv(0) * eU.getVCache().getDsf(0);
// upwinding
@@ -139,7 +137,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
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;
+ A1 += dmu * (rhoU - fluid.rho->eval(e, phi, k)) * fluid.mach->evalGrad(e, phi, k) * dSf.transpose() * dPhi * dPhi.transpose() * dSf * wdj;
}
}
return std::make_tuple(A1, A2);
@@ -155,7 +153,6 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
* + \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)
{
@@ -169,8 +166,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
{
// Density and gradients
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
- double rho = fluid.rho.eval(e, phi, k);
- double dRho = fluid.rho.evalD(e, phi, k);
+ double rho = fluid.rho->eval(e, phi, k);
+ double dRho = fluid.rho->evalGrad(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
Eigen::MatrixXd const &iJ = e.getJinv(k);
// Jacobian gradients
@@ -193,7 +190,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
}
// Supersonic contribution
Eigen::MatrixXd A2;
- double mach = fluid.mach.eval(e, phi, 0);
+ double mach = fluid.mach->eval(e, phi, 0);
if (mach > mCO)
{
// switching function and gradient
@@ -201,8 +198,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
double dmu = 2 * muC * mCO * mCO / (mach * mach * mach);
// upwind density and gradient
Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
- double rhoU = fluid.rho.eval(eU, phi, 0);
- double dRhoU = fluid.rho.evalD(eU, phi, 0);
+ double rhoU = fluid.rho->eval(eU, phi, 0);
+ double dRhoU = fluid.rho->evalGrad(eU, phi, 0);
Eigen::MatrixXd const &iJU = eU.getJinv(0);
std::vector<Eigen::MatrixXd> const &dJU = eU.getGradJ(0);
@@ -213,8 +210,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
{
// Density, Mach and Gradients
Eigen::VectorXd dPhi = e.computeGradient(phi, k);
- double rho = fluid.rho.eval(e, phi, k);
- double dM = fluid.mach.evalD(e, phi, k);
+ double rho = fluid.rho->eval(e, phi, k);
+ double dM = fluid.mach->evalGrad(e, phi, k);
Eigen::MatrixXd const &dSf = cache.getDsf(k);
Eigen::MatrixXd const &iJ = e.getJinv(k);
// Jacobian gradients
diff --git a/flow/src/wPotentialResidual.h b/flow/src/wPotentialResidual.h
index 62dde67b..79722718 100644
--- a/flow/src/wPotentialResidual.h
+++ b/flow/src/wPotentialResidual.h
@@ -27,6 +27,7 @@ namespace flow
/**
* @brief Formulation of nonlinear potential equation residuals
+ * @todo mu as F0ElMu?
* @todo split stabilization terms from subsonic terms?
*/
class FLOW_API PotentialResidual
diff --git a/flow/src/wProblem.cpp b/flow/src/wProblem.cpp
index 7d3d069f..46bc8ad6 100644
--- a/flow/src/wProblem.cpp
+++ b/flow/src/wProblem.cpp
@@ -29,96 +29,29 @@ using namespace flow;
Problem::Problem(std::shared_ptr<MshData> _msh, int dim, double aoa, double aos,
double minf, double sref, double cref, double xref,
double yref, double zref) : msh(_msh), nDim(dim), alpha(aoa),
- beta(aos), M_inf(minf), S_ref(sref),
- c_ref(cref)
+ beta(aos), M_inf(minf),
+ S_ref(sref), c_ref(cref),
+ x_ref(xref, yref, zref),
+ dirD(dim, sref, aoa, aos),
+ dirS(dim, sref, aoa, aos),
+ dirL(dim, sref, aoa, aos)
{
- x_ref(0) = xref;
- x_ref(1) = yref;
- 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)
- return Eigen::Vector2d(cos(alpha), sin(alpha));
- else
- 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));
+ obs.push_back(&dirD);
+ obs.push_back(&dirS);
+ obs.push_back(&dirL);
}
/**
* @brief Set the fluid medium
- * @authors Adrien Crovato
*/
void Problem::set(std::shared_ptr<Medium> m)
{
medium = m;
+ obs.push_back(m->phiInf);
}
/**
* @brief Add a boundary surface
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Boundary> b)
{
@@ -127,34 +60,33 @@ void Problem::add(std::shared_ptr<Boundary> b)
/**
* @brief Add initial condition
- * @authors Adrien Crovato
*/
void Problem::set(std::shared_ptr<Initial> i)
{
iIC = i;
+ obs.push_back(i->f);
}
/**
* @brief Add Dirichlet boundary condition
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Dirichlet> d)
{
dBCs.push_back(d);
+ obs.push_back(d->f);
}
/**
* @brief Add freestream (Neumann) boundary condition
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Freestream> f)
{
fBCs.push_back(f);
+ obs.push_back(f->f);
}
/**
* @brief Add wake boundary condition
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Wake> w)
{
@@ -163,7 +95,6 @@ void Problem::add(std::shared_ptr<Wake> w)
/**
* @brief Add Kutta condition
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Kutta> k)
{
@@ -172,13 +103,22 @@ void Problem::add(std::shared_ptr<Kutta> k)
/**
* @brief Add Blowing boundary condition
- * @authors Adrien Crovato
*/
void Problem::add(std::shared_ptr<Blowing> b)
{
bBCs.push_back(b);
}
+/**
+ * @brief Update the angle of attack
+ */
+void Problem::update(double aoa)
+{
+ alpha = aoa;
+ for (auto o : obs)
+ o->update(aoa);
+}
+
/**
* @brief Initialize the elements precomputed values
*/
@@ -186,7 +126,7 @@ void Problem::initElems()
{
// Update volume
for (auto e : medium->tag->elems)
- e->initValues(true);
+ e->initValues(true, 1);
// Update surface (boundaries)
for (auto surf : bnds)
for (auto e : surf->groups[0]->tag->elems)
@@ -212,7 +152,6 @@ void Problem::initElems()
/**
* @brief Initialize the elements precomputed gradients
- * @authors Adrien Crovato
*/
void Problem::initGradElems()
{
@@ -244,12 +183,11 @@ void Problem::initGradElems()
void Problem::pin()
{
if (dBCs.empty())
- this->add(std::make_shared<Dirichlet>(this->msh, this->fBCs.front()->tag->no, this->medium->phiInf, true));
+ this->add(std::make_shared<Dirichlet>(msh, fBCs.front()->tag->no, nDim, alpha, beta, true));
}
/**
* @brief Check that Problem is not empty and that elements are supported
- * @authors Adrien Crovato
*/
void Problem::check() const
{
diff --git a/flow/src/wProblem.h b/flow/src/wProblem.h
index 152aa1f2..53f2e5b3 100644
--- a/flow/src/wProblem.h
+++ b/flow/src/wProblem.h
@@ -19,6 +19,7 @@
#include "flow.h"
#include "wObject.h"
+#include "wF1Ct.h"
#include <memory>
#include <iostream>
#include <vector>
@@ -30,42 +31,45 @@ namespace flow
/**
* @brief Manage the problem
*
- * Contains mesh, medium, boundaries and reference value
+ * Contains freestream definition, reference values and physical groups
* @authors Adrien Crovato
*/
class FLOW_API Problem : public fwk::wSharedObject
{
public:
std::shared_ptr<tbox::MshData> msh; ///< Mesh structure
+ int nDim; ///< Problem dimension
+ // Freestream
+ double alpha; ///< Angle of attack
+ double beta; ///< Angle of sideslip
+ double M_inf; ///< Mach number
+ // Reference values
+ double S_ref; ///< Reference surface
+ double c_ref; ///< Reference chord
+ Eigen::Vector3d x_ref; ///< Reference center point (for moment computation)
#ifndef SWIG
- std::shared_ptr<Medium> medium; ///< Fluid
- std::vector<std::shared_ptr<Boundary>> bnds; ///< Boundaries
-
+ // Loads direction
+ F1CtDrag dirD; // drag
+ F1CtSide dirS; // sideforce
+ F1CtLift dirL; // lift
+ // Physical groups
+ std::shared_ptr<Medium> medium; ///< Fluid
+ std::vector<std::shared_ptr<Boundary>> bnds; ///< Boundaries
std::shared_ptr<Initial> iIC; ///< Initial condition
std::vector<std::shared_ptr<Dirichlet>> dBCs; ///< Dirichlet boundary conditions
std::vector<std::shared_ptr<Freestream>> fBCs; ///< Freestream boundary conditions
std::vector<std::shared_ptr<Wake>> wBCs; ///< Wake boundary condition
std::vector<std::shared_ptr<Kutta>> kSCs; ///< Kutta condition
std::vector<std::shared_ptr<Blowing>> bBCs; ///< Blowing (transpiration) boundary condition
+ // Observers
+ std::vector<fwk::Observer *> obs; ///< classes depending on problem variables
#endif
- // Reference values
- int nDim; ///< Problem dimension
- double alpha; ///< Angle of attack
- double beta; ///< Angle of sideslip
- double M_inf; ///< Freestream Mach number
- double S_ref; ///< Reference surface
- double c_ref; ///< Reference chord
- Eigen::Vector3d x_ref; ///< Reference center point (for moment computation)
Problem(std::shared_ptr<tbox::MshData> _msh, int dim, double aoa, double aos,
double minf, double sref, double cref, double xref, double yref,
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);
@@ -74,6 +78,7 @@ public:
void add(std::shared_ptr<Wake> w);
void add(std::shared_ptr<Kutta> k);
void add(std::shared_ptr<Blowing> b);
+ void update(double aoa);
#ifndef SWIG
void check() const;
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 7fae0493..991721ff 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -160,17 +160,11 @@ void Solver::save(MshExport *mshWriter, int n)
*/
void Solver::computeLoad()
{
- // Reset coefficients
+ // Reset coefficients and compute loads
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)
{
// Reset
@@ -180,7 +174,7 @@ void Solver::computeLoad()
for (auto e : sur->groups[0]->tag->elems)
{
// Build nodal pressure load
- Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, pbl->medium->cP);
+ 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)
@@ -201,9 +195,9 @@ void Solver::computeLoad()
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
- sur->Cd = Cf.dot(dirD);
- sur->Cs = Cf.dot(dirS);
- sur->Cl = Cf.dot(dirL);
+ sur->Cd = Cf.dot(pbl->dirD.eval());
+ sur->Cs = Cf.dot(pbl->dirS.eval());
+ sur->Cl = Cf.dot(pbl->dirL.eval());
// compute total
Cl += sur->Cl;
Cd += sur->Cd;
@@ -223,7 +217,7 @@ void Solver::computeFlow()
tbb::parallel_do(fluid->neMap.begin(), fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
Node *n = p.first;
// Perturbation potential
- vPhi[n->row] = phi[n->row] - fluid->phiInf.eval(n->pos);
+ vPhi[n->row] = phi[n->row] - fluid->phiInf->eval(n->pos);
// Velocity (averaged at nodes)
U[n->row] = Eigen::Vector3d::Zero();
for (auto e : p.second)
@@ -236,17 +230,17 @@ void Solver::computeFlow()
// Density (averaged at nodes)
rho[n->row] = 0.;
for (auto e : p.second)
- rho[n->row] += fluid->rho.eval(*e, phi, 0);
+ rho[n->row] += fluid->rho->eval(*e, phi, 0);
rho[n->row] /= p.second.size();
// Mach number (averaged at nodes)
M[n->row] = 0.;
for (auto e : p.second)
- M[n->row] += fluid->mach.eval(*e, phi, 0);
+ M[n->row] += fluid->mach->eval(*e, phi, 0);
M[n->row] /= p.second.size();
// Pressure coefficient (averaged at nodes)
Cp[n->row] = 0.;
for (auto e : p.second)
- Cp[n->row] += fluid->cP.eval(*e, phi, 0);
+ Cp[n->row] += fluid->cP->eval(*e, phi, 0);
Cp[n->row] /= p.second.size();
});
diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp
index 12562c7c..fcfd06b7 100644
--- a/flow/src/wWakeResidual.cpp
+++ b/flow/src/wWakeResidual.cpp
@@ -30,8 +30,9 @@ using namespace flow;
* 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
+ * NB: v_inf = [1, 0, 0]
*/
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElement const &we, std::vector<double> const &phi)
{
// Get shape functions and Gauss points
Cache &sCache = we.surUpE->getVCache();
@@ -51,7 +52,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
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(we.surUpE->getVol()) * dSf * vUpJ.transpose() * vi;
+ dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * vUpJ.transpose() * Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1);
// Velocity
Eigen::RowVectorXd Aup = we.volUpE->computeGradient(phi, 0).transpose() * vUpJ * vUpDsf;
Eigen::RowVectorXd Alw = we.volLwE->computeGradient(phi, 0).transpose() * vLwJ * vLwDsf;
@@ -78,8 +79,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
*
* 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
+ * NB: v_inf = [1, 0, 0]
*/
-std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement const &we, std::vector<double> const &phi)
{
// Get shape functions and Gauss points
Cache &sCache = we.surUpE->getVCache();
@@ -95,7 +97,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
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(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * vi;
+ dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1);
// Velocity squared
Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
@@ -124,12 +126,13 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
*
* 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
+ * NB: v_inf = [1, 0, 0]
*/
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(WakeElement const &we, std::vector<double> const &phi)
{
// Build
Eigen::MatrixXd Au, Al;
- std::tie(Au, Al) = WakeResidual::buildFixed(we, phi, vi);
+ std::tie(Au, Al) = WakeResidual::buildFixed(we, phi);
return std::make_tuple(2 * Au, 2 * Al);
}
@@ -140,8 +143,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientFlow(Wak
* = \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
+ * NB: v_inf = [1, 0, 0]
*/
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildGradientMesh(WakeElement const &we, std::vector<double> const &phi)
{
// Get shape functions and Gauss points
Cache &sCache = we.surUpE->getVCache();
@@ -162,6 +166,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
size_t nSur = we.surUpE->nodes.size();
// Stabilization term
+ Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).block(0, 0, nDim, 1);
double sqSurf = sqrt(we.surUpE->getVol());
Eigen::MatrixXd dSf(nRow, nDim);
for (size_t i = 0; i < nRow; ++i)
diff --git a/flow/src/wWakeResidual.h b/flow/src/wWakeResidual.h
index fc536849..5bad9298 100644
--- a/flow/src/wWakeResidual.h
+++ b/flow/src/wWakeResidual.h
@@ -32,12 +32,12 @@ 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);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildFixed(WakeElement const &we, std::vector<double> const &phi);
// 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);
+ static std::tuple<Eigen::VectorXd, Eigen::VectorXd> build(WakeElement const &we, std::vector<double> const &phi);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(WakeElement const &we, std::vector<double> const &phi);
// Adjoint
- static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(WakeElement const &we, std::vector<double> const &phi, Eigen::VectorXd const &vi);
+ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(WakeElement const &we, std::vector<double> const &phi);
};
} // namespace flow
diff --git a/flow/tests/adjoint.py b/flow/tests/adjoint.py
index e0134b0f..990fb526 100644
--- a/flow/tests/adjoint.py
+++ b/flow/tests/adjoint.py
@@ -43,7 +43,6 @@ def main():
# define flow variables
alpha = 2*np.pi/180
M_inf = 0.7
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
@@ -58,7 +57,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/adjoint3.py b/flow/tests/adjoint3.py
index 670ed6bf..c52cf8d7 100644
--- a/flow/tests/adjoint3.py
+++ b/flow/tests/adjoint3.py
@@ -43,7 +43,6 @@ def main():
# define flow variables
alpha = 3*np.pi/180
M_inf = 0.3
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
span = 1
dim = 3
@@ -59,7 +58,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/bli.py b/flow/tests/bli.py
index 73728691..3c7b3584 100644
--- a/flow/tests/bli.py
+++ b/flow/tests/bli.py
@@ -55,7 +55,6 @@ def main():
alpha = 5*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.5
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
tol = 1e-4 # tolerance for the VII (usually one drag count)
@@ -69,7 +68,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
tms['pre'].stop()
# solve viscous problem
diff --git a/flow/tests/cylinder.py b/flow/tests/cylinder.py
index 2a0b51c4..1cd3a0d9 100644
--- a/flow/tests/cylinder.py
+++ b/flow/tests/cylinder.py
@@ -40,7 +40,6 @@ def main():
alpha = 3*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.1
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 2
# define dimension and mesh size
@@ -60,7 +59,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True)
+ pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True)
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/cylinder2D5.py b/flow/tests/cylinder2D5.py
index 14f943df..3999b7e4 100644
--- a/flow/tests/cylinder2D5.py
+++ b/flow/tests/cylinder2D5.py
@@ -40,7 +40,6 @@ def main():
alpha = 3*math.pi/180
U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
M_inf = 0.1
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension and mesh size
@@ -62,7 +61,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True)
+ pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True)
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/cylinder3.py b/flow/tests/cylinder3.py
index 2a4b36c5..fe8328b0 100644
--- a/flow/tests/cylinder3.py
+++ b/flow/tests/cylinder3.py
@@ -41,7 +41,6 @@ def main():
alpha = 3*math.pi/180
U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
M_inf = 0.1
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension and mesh size
@@ -65,7 +64,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True)
+ pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True)
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/lift.py b/flow/tests/lift.py
index a5757ec9..52cce3bf 100644
--- a/flow/tests/lift.py
+++ b/flow/tests/lift.py
@@ -41,7 +41,6 @@ def main():
# define flow variables
alpha = 2*math.pi/180
M_inf = 0.7
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
@@ -54,7 +53,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/lift3.py b/flow/tests/lift3.py
index da872d0f..cd05acaf 100644
--- a/flow/tests/lift3.py
+++ b/flow/tests/lift3.py
@@ -40,7 +40,6 @@ def main():
# define flow variables
alpha = 3*math.pi/180
M_inf = 0.3
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension and mesh size
@@ -62,7 +61,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
diff --git a/flow/tests/meshDef.py b/flow/tests/meshDef.py
index d258e5a6..f6fdba13 100644
--- a/flow/tests/meshDef.py
+++ b/flow/tests/meshDef.py
@@ -41,7 +41,6 @@ def main():
# define flow variables
alpha0 = 0 # must be zero for this testfile
M_inf = 0.3
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
# parameters for mesh deformation
@@ -64,9 +63,9 @@ 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+dx, dy, 0., 'airfoil', te = 'te')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, 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')
+ pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te')
tms['pre'].stop()
# solver problem for 0°AOA
diff --git a/flow/tests/meshDef3.py b/flow/tests/meshDef3.py
index b1e41191..11459bb6 100644
--- a/flow/tests/meshDef3.py
+++ b/flow/tests/meshDef3.py
@@ -40,7 +40,6 @@ def main():
# define flow variables
alpha0 = 0 # must be zero for this testfile
M_inf = 0.3
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension and mesh size
@@ -69,7 +68,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, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solver problem for 0°AOA
diff --git a/flow/tests/nonlift.py b/flow/tests/nonlift.py
index ef6c22cc..15e2bbf2 100644
--- a/flow/tests/nonlift.py
+++ b/flow/tests/nonlift.py
@@ -42,7 +42,6 @@ def main():
# define flow variables
alpha = 0*math.pi/180 # must be zero for this testfile
M_inf = 0.8
- M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
@@ -55,7 +54,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
tms['pre'].stop()
# solve problem
diff --git a/flow/validation/agard.py b/flow/validation/agard.py
index edc13ba3..8df7bea8 100644
--- a/flow/validation/agard.py
+++ b/flow/validation/agard.py
@@ -70,7 +70,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
@@ -100,6 +100,8 @@ def main():
tests = CTests()
tests.add(CTest('CL', solver.Cl, 0.062, 1e-1))
tests.add(CTest('CD', solver.Cd, 0.0006, 1e-1))
+ tests.add(CTest('CS', solver.Cs, 0.0017, 1e-1))
+ tests.add(CTest('CM', solver.Cm, -0.064, 1e-1))
tests.run()
# eof
diff --git a/flow/validation/onera.py b/flow/validation/onera.py
index ab71b2ba..8d742f02 100644
--- a/flow/validation/onera.py
+++ b/flow/validation/onera.py
@@ -43,7 +43,6 @@ def main():
# define flow variables
alpha = 3.06*np.pi/180
M_inf = 0.839
- M_crit = 5 # Squared critical Mach number (above which density is modified)
dim = 3
# define dimension and mesh size
@@ -68,7 +67,7 @@ def main():
# set the problem and add medium, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
diff --git a/fwk/src/fwk.h b/fwk/src/fwk.h
index dff0dee4..85c99132 100644
--- a/fwk/src/fwk.h
+++ b/fwk/src/fwk.h
@@ -48,6 +48,7 @@ namespace fwk
class wObject;
class Extract;
+class Observer;
class HMSTime;
class Time;
diff --git a/fwk/src/wObserver.cpp b/fwk/src/wObserver.cpp
new file mode 100644
index 00000000..0162b4a3
--- /dev/null
+++ b/fwk/src/wObserver.cpp
@@ -0,0 +1,25 @@
+/*
+ * 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 "wObserver.h"
+#include <stdexcept>
+
+using namespace fwk;
+
+void Observer::update(double data)
+{
+ throw std::runtime_error("fwk::Observer::update() not implemented!\n");
+}
\ No newline at end of file
diff --git a/fwk/src/wObserver.h b/fwk/src/wObserver.h
new file mode 100644
index 00000000..cabc058c
--- /dev/null
+++ b/fwk/src/wObserver.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 WOBSERVER_H
+#define WOBSERVER_H
+
+#include "fwk.h"
+
+namespace fwk
+{
+
+/**
+ * @brief Base class that obsverves another class, and that updates accordingly
+ */
+
+class FWK_API Observer
+{
+public:
+ Observer() {}
+ virtual ~Observer() {}
+ virtual void update(double data);
+};
+
+} // namespace fwk
+
+#endif //WOBSERVER_H
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 35bab22c..7bca3652 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -71,7 +71,6 @@ Solver::~Solver()
}
/**
- * @todo detecter/traiter les jacos negatifs pour tous les éléments
* @todo planter clairement si pas de Medium
* @todo calculer la taille du systeme en fct du nombre de noeuds des physical groups utilisés
*/
diff --git a/tbox/src/wFct0.cpp b/tbox/src/wFct0.cpp
index b0a7b0dc..2f4618d3 100644
--- a/tbox/src/wFct0.cpp
+++ b/tbox/src/wFct0.cpp
@@ -66,13 +66,6 @@ void PwLf::write(std::ostream &out) 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");
-}
-
-//-------------------------------------------------------------------------
-
void Fct0C::write(std::ostream &out) const
{
out << "Fct0C=" << v;
diff --git a/tbox/src/wFct0.h b/tbox/src/wFct0.h
index 81da6cd4..eb6d4166 100644
--- a/tbox/src/wFct0.h
+++ b/tbox/src/wFct0.h
@@ -41,10 +41,9 @@ public:
/**
* @brief Scalar functions to be integrated over an element
- *
+ *
* eval() evaluates the function at kth Gauss point
- * evalD() evaluates the derivative at the kth Gauss point
- * @authors Romain Boman, Adrien Crovato
+ * @authors Romain Boman
*/
class TBOX_API Fct0 : public fwk::wObject
{
@@ -62,9 +61,6 @@ public:
virtual double eval(Element const &e,
std::vector<double> const &u,
size_t k) const = 0;
- virtual double evalD(Element const &e,
- std::vector<double> const &u,
- size_t k) const;
#endif
};
@@ -78,7 +74,6 @@ class TBOX_API Fct0C : public Fct0
public:
Fct0C(double _v) : Fct0(), v(_v) {}
- void update(double _v) { v = _v; }
#ifndef SWIG
virtual double eval(Element const &e,
@@ -87,12 +82,6 @@ public:
{
return v;
}
- virtual double evalD(Element const &e,
- std::vector<double> const &u,
- size_t k) const override
- {
- return 0;
- }
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/tbox/src/wFct1.h b/tbox/src/wFct1.h
index 439747f2..0a3c546d 100644
--- a/tbox/src/wFct1.h
+++ b/tbox/src/wFct1.h
@@ -27,7 +27,7 @@ namespace tbox
/**
* @brief Vector functions to be integrated over an element
- *
+ *
* eval() evaluates the function at kth Gauss point
* @authors Romain Boman, Adrien Crovato
*/
@@ -47,9 +47,6 @@ public:
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
};
@@ -76,12 +73,6 @@ public:
{
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
};
--
GitLab