From 685c578c710741d4a9e92b38a59d3421b03a259e Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Sun, 19 Dec 2021 15:52:12 +0100
Subject: [PATCH] Add Kutta for 3D computations. Adapt geos, API and
computations. Refactor tests. Minor bugfix.
---
dart/api/core.py | 11 +-
dart/api/internal/trim.py | 2 +-
dart/benchmark/onera.py | 2 +-
dart/benchmark/onera_lfs.msh | 4 +-
dart/cases/coyote.py | 5 +-
dart/cases/n0012_3.py | 73 -----
dart/cases/wbht.py | 7 +-
dart/cases/wht.py | 4 +-
dart/default.py | 17 +-
dart/models/agard445.geo | 2 +-
dart/models/coyote.geo | 6 +-
dart/models/cylinder_2D5.geo | 161 ----------
dart/models/cylinder_3.geo | 203 ------------
dart/models/n0012_3.geo | 442 --------------------------
dart/models/oneraM6.geo | 2 +-
dart/models/rae_25.geo | 400 +++++++++++++++++++++++
dart/models/rae_3.geo | 430 +++++++++++++++++++++++++
dart/models/wbht.geo | 12 +-
dart/models/wht.geo | 10 +-
dart/src/wAdjoint.cpp | 4 +-
dart/src/wKuttaElement.cpp | 37 ++-
dart/src/wKuttaElement.h | 21 +-
dart/src/wKuttaResidual.cpp | 101 ++++--
dart/src/wProblem.cpp | 99 ++----
dart/src/wProblem.h | 4 +
dart/src/wSolver.cpp | 3 +
dart/src/wWakeElement.h | 4 +-
dart/src/wWakeResidual.cpp | 10 +-
dart/tests/adjoint.py | 16 +-
dart/tests/bli.py | 4 +-
dart/tests/cylinder.py | 2 +-
dart/tests/cylinder2D5.py | 134 --------
dart/tests/cylinder3.py | 137 --------
dart/tests/lift.py | 3 +-
dart/tests/lift3.py | 101 ------
dart/tests/meshDef3.py | 134 --------
dart/tests/{meshDef.py => morpher.py} | 10 +-
dart/tests/nonlift.py | 3 +-
dart/tests/{adjoint3.py => rae_25.py} | 80 ++---
dart/tests/rae_3.py | 161 ++++++++++
dart/validation/agard.py | 2 +-
dart/validation/onera.py | 2 +-
ext/amfe | 2 +-
43 files changed, 1277 insertions(+), 1590 deletions(-)
delete mode 100644 dart/cases/n0012_3.py
delete mode 100644 dart/models/cylinder_2D5.geo
delete mode 100644 dart/models/cylinder_3.geo
delete mode 100644 dart/models/n0012_3.geo
create mode 100644 dart/models/rae_25.geo
create mode 100644 dart/models/rae_3.geo
delete mode 100644 dart/tests/cylinder2D5.py
delete mode 100644 dart/tests/cylinder3.py
delete mode 100644 dart/tests/lift3.py
delete mode 100644 dart/tests/meshDef3.py
rename dart/tests/{meshDef.py => morpher.py} (96%)
rename dart/tests/{adjoint3.py => rae_25.py} (62%)
create mode 100644 dart/tests/rae_3.py
diff --git a/dart/api/core.py b/dart/api/core.py
index 7da8be7..e729b8d 100644
--- a/dart/api/core.py
+++ b/dart/api/core.py
@@ -136,7 +136,7 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
if 'Symmetry' in params:
mshCrck.addBoundaries([params['Symmetry']])
if 'WakeTips' in params:
- mshCrck.setExcluded(params['WakeTips'][i]) # 2.5D computations
+ mshCrck.setExcluded(params['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
mshCrck.run()
# save the updated mesh in native (gmsh) format and then set the writer to the requested format
tbox.GmshExport(_msh).save(_msh.name)
@@ -191,12 +191,19 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']])
_pbl.add(bnd)
# add Wake/Kutta boundary conditions
+ # TODO refactor Wake "exclusion" for 3D?
if _dim == 2:
_pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
_pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
else:
for i in range(len(params['Wakes'])):
- _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['TeTips'][i]]))
+ if 'WakeExs' in params:
+ _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeExs'][i]])) # 3D + fuselage
+ elif 'WakeTips' in params:
+ _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeTips'][i]])) # 3D
+ else:
+ _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid']])) # 2.5D
+ _pbl.add(dart.Kutta(_msh, [params['Tes'][i], params['Wakes'][i]+'_', params['Wings'][i], params['Fluid']]))
# Direct (forward) solver creation
# NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index 2555218..be1c330 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -111,7 +111,7 @@ class Trim:
self.sol.save(self.wrtr)
# extract Cp
if self.dim == 2:
- Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
+ Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "")
elif self.dim == 3 and self.format == 'vtk' and self.slice:
dU.writeSlices(self.msh.name, self.slice, self.tag)
diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py
index 08e2f85..75209e2 100644
--- a/dart/benchmark/onera.py
+++ b/dart/benchmark/onera.py
@@ -66,7 +66,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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 = 'wakeTip')
tms['pre'].stop()
# solve problem
diff --git a/dart/benchmark/onera_lfs.msh b/dart/benchmark/onera_lfs.msh
index bad8834..b69a997 100644
--- a/dart/benchmark/onera_lfs.msh
+++ b/dart/benchmark/onera_lfs.msh
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:da782b51f3be6ff2b8d222ed9ad1aba46849b2d018df294ac5f780a0258a74c8
-size 39450759
+oid sha256:c235746b196f5b428da900f56a024eb69c5d1baa2766409c0e1413109c2cd7fb
+size 36951138
diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py
index 6dda2e0..fc493e3 100644
--- a/dart/cases/coyote.py
+++ b/dart/cases/coyote.py
@@ -44,7 +44,7 @@ def getParam():
p['Dim'] = 3 # Problem dimension
p['Format'] = 'vtk' # Save format (vtk or gmsh)
p['Slice'] = [1.8, 3.5, 5.3] # array of (y) coordinates to perform slice along the span (empty if none)
- p['TagId'] = 9 # tag number of physical group to be sliced
+ p['TagId'] = 11 # tag number of physical group to be sliced
# Markers
p['Fluid'] = 'field' # Name of physical group containing the fluid
p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
@@ -52,7 +52,8 @@ def getParam():
p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
- p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
+ p['WakeExs'] = ['wakeExW', 'wakeExT'] # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
+ p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge
# Freestream
p['M_inf'] = 0.8 # Freestream Mach number
p['AoA_begin'] = 1 # Freestream angle of attack (begin)
diff --git a/dart/cases/n0012_3.py b/dart/cases/n0012_3.py
deleted file mode 100644
index 57700c2..0000000
--- a/dart/cases/n0012_3.py
+++ /dev/null
@@ -1,73 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## NACA 0012 wing configuration file for flow polar script
-# Adrien Crovato
-
-def getParam():
- import os.path
- from fwk.wutils import parseargs
- p = {}
- # Specific
- span = 1. # span length
- # Arguments
- args = parseargs()
- p['Threads'] = args.k
- p['Verb'] = args.verb
- # Input/Output
- p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012_3.geo') # Input file containing the model
- p['Pars'] = {'spn' : span, 'lgt' : 3., 'wdt' : 3., 'hgt' : 3., 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} # Parameters for input file model
- p['Dim'] = 3 # Problem dimension
- p['Format'] = 'vtk' # Save format (vtk or gmsh)
- p['Slice'] = [0.01, 0.25, 0.5, 0.75, 0.95] # array of (y) coordinates to perform slice along the span (empty if none)
- p['TagId'] = 4 # tag number of physical group to be sliced
- # Markers
- p['Fluid'] = 'field' # Name of physical group containing the fluid
- p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
- p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
- p['Wings'] = ['wing'] # LIST of names of physical groups containing the lifting surface body
- p['Wakes'] = ['wake'] # LIST of names of physical group containing the wake
- p['WakeTips'] = ['wakeTip'] # LIST of names of physical group containing the edge of the wake
- p['TeTips'] = ['teTip'] # LIST of names of physical group containing the edge of the wake and the trailing edge
- # Freestream
- p['M_inf'] = 0.3 # Freestream Mach number
- p['AoA_begin'] = -3. # Freestream angle of attack (begin)
- p['AoA_end'] = 3. # Freestream angle of attack (end)
- p['AoA_step'] = 1. # Freestream angle of attack (step)
- # Geometry
- p['S_ref'] = 1.*span # Reference surface length
- p['c_ref'] = 1. # Reference chord length
- p['x_ref'] = 0.25 # Reference point for moment computation (x)
- p['y_ref'] = 0. # Reference point for moment computation (y)
- p['z_ref'] = 0. # Reference point for moment computation (z)
- # Numerical
- p['LSolver'] = 'MUMPS' # Linear solver (PARDISO, GMRES, or MUMPS)
- p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
- p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
- p['Max_it'] = 10 # Solver maximum number of iterations
- return p
-
-def main():
- # Script call
- import dart.api.internal.polar as p
- polar = p.Polar(getParam())
- polar.run()
- polar.disp()
-
-if __name__ == "__main__":
- main()
diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py
index 4744d14..73a8f1d 100644
--- a/dart/cases/wbht.py
+++ b/dart/cases/wbht.py
@@ -56,7 +56,7 @@ def getParam():
p['Dim'] = 3 # Problem dimension
p['Format'] = 'vtk' # Save format (vtk or gmsh)
p['Slice'] = [0.83, 3.0, 4.0, 6.2, 6.9] # array of (y) coordinates to perform slice along the span (empty if none)
- p['TagId'] = 10 # tag number of physical group to be sliced
+ p['TagId'] = 12 # tag number of physical group to be sliced
# Markers
p['Fluid'] = 'field' # Name of physical group containing the fluid
p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
@@ -65,9 +65,10 @@ def getParam():
p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
- p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
+ p['WakeExs'] = ['wakeExW', 'wakeExT'] # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
+ p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge
# Freestream
- p['M_inf'] = 0.75 # Freestream Mach number
+ p['M_inf'] = 0.72 # Freestream Mach number
p['CL'] = 0.4 # Target lift coefficient
p['AoA'] = 0. # Freestream angle of attack (guess)
p['dCL'] = 7.5 # Slope of CL w.r.t AoA [1/rad] (guess)
diff --git a/dart/cases/wht.py b/dart/cases/wht.py
index ce73017..dd0fc7d 100644
--- a/dart/cases/wht.py
+++ b/dart/cases/wht.py
@@ -59,9 +59,9 @@ def getParam():
p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
- p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
+ p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge
# Freestream
- p['M_inf'] = 0.75 # Freestream Mach number
+ p['M_inf'] = 0.72 # Freestream Mach number
p['CL'] = 0.4 # Target lift coefficient
p['AoA'] = 0. # Freestream angle of attack (guess)
p['dCL'] = 7.5 # Slope of CL w.r.t AoA [1/rad] (guess)
diff --git a/dart/default.py b/dart/default.py
index 3a03df5..93cb86c 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -57,7 +57,7 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
gmshWriter.save(msh.name)
return msh, gmshWriter
-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):
+def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = 'te', tp = None, vsc = False, dbc = False):
"""Initialize problem
Parameters
@@ -97,7 +97,7 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld
te : str (optional)
name of physical tag contaning the trailing edge
tp : str (optional)
- name of physical tag contaning the wake tip (and the trailing edge for 3D)
+ name of physical tag contaning the wake tip (only for 3D)
vsc : str (optional)
name of physical tag contaning the body to apply transpiration B.C.
dbc : bool (optional)
@@ -121,15 +121,12 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld
pbl.add(dart.Freestream(msh, ff, dim, alpha, beta))
pbl.add(dart.Freestream(msh, dnstr, dim, alpha, beta))
# add wake and Kutta conditions
- if dim == 2 and te is not None:
- wake = dart.Wake(msh, [wk, wk+'_', fld])
- pbl.add(wake)
- pbl.add(dart.Kutta(msh, [te, wk+'_', sur, fld]))
- elif dim == 3:
- wake = dart.Wake(msh, [wk, wk+'_', fld, tp])
- pbl.add(wake)
+ if tp is None:
+ wake = dart.Wake(msh, [wk, wk+'_', fld]) # 2 and 2.5D
else:
- wake = None
+ wake = dart.Wake(msh, [wk, wk+'_', fld, tp]) # 3D
+ pbl.add(wake)
+ pbl.add(dart.Kutta(msh, [te, wk+'_', sur, fld]))
# add body of interest
bnd = dart.Body(msh, [sur, fld])
pbl.add(bnd)
diff --git a/dart/models/agard445.geo b/dart/models/agard445.geo
index c700751..d268664 100644
--- a/dart/models/agard445.geo
+++ b/dart/models/agard445.geo
@@ -354,7 +354,7 @@ Mesh.SmoothNormals = 1;
// Trailing edge and wake tip
Physical Line("wakeTip") = {162};
-Physical Line("teTip") = {61, 162};
+Physical Line("te") = {61};
// Internal Field:
Physical Volume("field") = {1};
diff --git a/dart/models/coyote.geo b/dart/models/coyote.geo
index c9e3615..bdf9a0c 100644
--- a/dart/models/coyote.geo
+++ b/dart/models/coyote.geo
@@ -61,9 +61,11 @@ Coherence;
/* Physical */
// Tips
-Physical Line("teTipW") = {155:158, 160:161, 164:167, 169:170};
+Physical Line("wakeExW") = {156:158, 160, 165:167, 169};
+Physical Line("teW") = {155, 161, 164, 170};
Physical Line("wakeTipW") = {160, 169};
-Physical Line("teTipT") = {173:176, 178:181};
+Physical Line("wakeExT") = {173, 175, 176, 178, 180, 181};
+Physical Line("teT") = {174, 179};
Physical Line("wakeTipT") = {173, 178};
// Downstream
diff --git a/dart/models/cylinder_2D5.geo b/dart/models/cylinder_2D5.geo
deleted file mode 100644
index cc70813..0000000
--- a/dart/models/cylinder_2D5.geo
+++ /dev/null
@@ -1,161 +0,0 @@
-/* 2D5 Cylinder */
-
-// Parameters
-// domain
-DefineConstant[ dmt = { 1.0, Name "Cylinder diameter" } ];
-DefineConstant[ lgt = { 3.0, Name "Channel length" } ];
-DefineConstant[ wdt = { 2.0, Name "Channel width" } ];
-DefineConstant[ hgt = { 3.0, Name "Channel height" } ];
-DefineConstant[ msN = { 0.1, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" } ];
-DefineConstant[ msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ];
-
-// mesh algo
-msh3D = 2;
-msh2D = 5;
-
-// Geometry
-// cylinder
-Point(0) = {0, 0, 0};
-Point(1) = {dmt/2, 0, 0, msN};
-Point(2) = {0, 0, dmt/2, msN};
-Point(3) = {-dmt/2, 0, 0, msN};
-Point(4) = {0, 0, -dmt/2, msN};
-Point(10) = {0, wdt, 0};
-Point(11) = {dmt/2, wdt, 0, msN};
-Point(12) = {0, wdt, dmt/2, msN};
-Point(13) = {-dmt/2, wdt, 0, msN};
-Point(14) = {0, wdt, -dmt/2, msN};
-
-// box
-Point(31) = {lgt/2, 0, hgt/2, msF};
-Point(32) = {-lgt/2, 0, hgt/2, msF};
-Point(33) = {-lgt/2, 0, -hgt/2, msF};
-Point(34) = {lgt/2, 0, -hgt/2, msF};
-Point(35) = {lgt/2, wdt, hgt/2, msF};
-Point(36) = {-lgt/2, wdt, hgt/2, msF};
-Point(37) = {-lgt/2, wdt, -hgt/2, msF};
-Point(38) = {lgt/2, wdt, -hgt/2, msF};
-
-// midplane
-Point(41) = {lgt/2, 0, 0, msF};
-Point(42) = {lgt/2, wdt, 0, msF};
-Point(43) = {-lgt/2, wdt, 0, msF};
-Point(44) = {-lgt/2, 0, 0, msF};
-
-// wing
-Circle(1) = {1, 0, 2};
-Circle(2) = {2, 0, 3};
-Circle(3) = {3, 0, 4};
-Circle(4) = {4, 0, 1};
-Circle(5) = {11, 10, 12};
-Circle(6) = {12, 10, 13};
-Circle(7) = {13, 10, 14};
-Circle(8) = {14, 10, 11};
-Line(9) = {1, 11};
-Line(10) = {2, 12};
-Line(11) = {3, 13};
-Line(12) = {4, 14};
-
-// box
-Line(21) = {41, 31};
-Line(22) = {31, 32};
-Line(23) = {32, 44};
-Line(24) = {44, 33};
-Line(25) = {33, 34};
-Line(26) = {34, 41};
-Line(27) = {42, 35};
-Line(28) = {35, 36};
-Line(29) = {36, 43};
-Line(30) = {43, 37};
-Line(31) = {37, 38};
-Line(32) = {38, 42};
-Line(33) = {31, 35};
-Line(34) = {32, 36};
-Line(35) = {33, 37};
-Line(36) = {34, 38};
-
-// midplane
-Line(41) = {1, 41};
-Line(42) = {11, 42};
-Line(43) = {13, 43};
-Line(44) = {3, 44};
-Line(45) = {41, 42};
-Line(46) = {42, 11};
-Line(47) = {13, 43};
-Line(48) = {43, 44};
-
-// wing
-Line Loop(1) = {1, 10, -5, -9};
-Surface(1) = {-1};
-Line Loop(2) = {2, 11, -6, -10};
-Surface(2) = {-2};
-Line Loop(3) = {3, 12, -7, -11};
-Surface(3) = {-3};
-Line Loop(4) = {4, 9, -8, -12};
-Surface(4) = {-4};
-
-// symmetry
-Line Loop(5) = {41, 21, 22, 23, -44, -2, -1};
-Plane Surface(5) = {-5};
-Line Loop(12) = {42, 27, 28, 29, -43, -6, -5};
-Plane Surface(12) = {12};
-Line Loop(6) = {44, 24, 25, 26, -41, -4, -3};
-Plane Surface(6) = {-6};
-Line Loop(13) = {43, 30, 31, 32, -42, -8, -7};
-Plane Surface(13) = {13};
-
-// upstream
-Line Loop(7) = {48, -23, 34, 29};
-Plane Surface(7) = {-7};
-Line Loop(8) = {48, 24, 35, -30};
-Plane Surface(8) = {8};
-
-// downstream
-Line Loop(9) = {45, 27, -33, -21};
-Plane Surface(9) = {-9};
-Line Loop(10) = {45, -32, -36, 26};
-Plane Surface(10) = {10};
-
-// farfield
-Line Loop(11) = {22, 34, -28, -33};
-Plane Surface(11) = {11};
-Line Loop(14) = {25, 36, -31, -35};
-Plane Surface(14) = {14};
-
-// wake
-Line Loop(15) = {41, 45, -42, -9};
-Plane Surface(15) = {15};
-
-// internal
-Line Loop(16) = {11, 43, 48, -44};
-Plane Surface(16) = {16};
-
-// field
-Surface Loop(1) = {7, 5, 9, 12, 11, 2, 1, 16, 15};
-Volume(1) = {1};
-Surface Loop(2) = {6, 8, 14, 10, 13, 4, 3, 16, 15};
-Volume(2) = {2};
-
-// Mesh
-// 2D
-MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} = msh2D;
-// 3D
-Mesh.Algorithm3D = msh3D;
-Mesh.OptimizeNetgen = 1;
-Mesh.Smoothing = 10;
-Mesh.SmoothNormals = 1;
-
-// Physical
-Physical Line("te") = {9};
-Physical Surface("cylinder") = {1, 2};
-Physical Surface("cylinder_") = {3, 4};
-Physical Surface("symmetry") = {5, 12};
-Physical Surface("symmetry_") = {6, 13};
-Physical Surface("upstream") = {7, 8};
-Physical Surface("downstream") = {9};
-Physical Surface("downstream_") = {10};
-Physical Surface("farfield") = {11, 14};
-Physical Surface("wake") = {15};
-Physical Volume("field") = {1};
-Physical Volume("field_") = {2};
-
diff --git a/dart/models/cylinder_3.geo b/dart/models/cylinder_3.geo
deleted file mode 100644
index f5c9602..0000000
--- a/dart/models/cylinder_3.geo
+++ /dev/null
@@ -1,203 +0,0 @@
-/* 3D Cylinder */
-
-// Parameters
-// domain
-DefineConstant[ dmt = { 1.0, Name "Cylinder diameter" } ];
-DefineConstant[ spn = { 2.5, Name "Cylinder span" } ];
-DefineConstant[ tpt = { .05, Name "Cylinder tip thickness" } ]; // percentage of span
-DefineConstant[ lgt = { 8.0, Name "Channel length" } ];
-DefineConstant[ wdt = { 4.0, Name "Channel width" } ];
-DefineConstant[ hgt = { 6.0, Name "Channel height" } ];
-DefineConstant[ msN = { 0.1, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" } ];
-DefineConstant[ msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ];
-
-// mesh algo
-msh3D = 2;
-msh2D = 5;
-
-// Geometry
-// cylinder
-Point(0) = {0, 0, 0};
-Point(1) = {dmt/2, 0, 0, msN};
-Point(2) = {0, 0, dmt/2, msN};
-Point(3) = {-dmt/2, 0, 0, msN};
-Point(4) = {0, 0, -dmt/2, msN};
-Point(10) = {0, spn, 0};
-Point(11) = {dmt/2, spn, 0, msN};
-Point(12) = {0, spn, dmt/2, msN};
-Point(13) = {-dmt/2, spn, 0, msN};
-Point(14) = {0, spn, -dmt/2, msN};
-
-// tip
-Point(21) = {0, spn+tpt*spn, 0, msN};
-
-// box
-Point(31) = {lgt/2, 0, hgt/2, msF};
-Point(32) = {-lgt/2, 0, hgt/2, msF};
-Point(33) = {-lgt/2, 0, -hgt/2, msF};
-Point(34) = {lgt/2, 0, -hgt/2, msF};
-Point(35) = {lgt/2, wdt, hgt/2, msF};
-Point(36) = {-lgt/2, wdt, hgt/2, msF};
-Point(37) = {-lgt/2, wdt, -hgt/2, msF};
-Point(38) = {lgt/2, wdt, -hgt/2, msF};
-
-// midplane
-Point(41) = {lgt/2, 0, 0, msF};
-Point(42) = {lgt/2, spn, 0, msF};
-Point(43) = {lgt/2, wdt, 0, msF};
-Point(44) = {dmt/2, wdt, 0, msF};
-Point(45) = {0, wdt, 0, msF};
-Point(46) = {-dmt/2, wdt, 0, msF};
-Point(47) = {-lgt/2, wdt, 0, msF};
-Point(48) = {-lgt/2, spn, 0, msF};
-Point(49) = {-lgt/2, 0, 0, msF};
-
-// wing
-Circle(1) = {1, 0, 2};
-Circle(2) = {2, 0, 3};
-Circle(3) = {3, 0, 4};
-Circle(4) = {4, 0, 1};
-Circle(11) = {11, 10, 12};
-Circle(12) = {12, 10, 13};
-Circle(13) = {13, 10, 14};
-Circle(14) = {14, 10, 11};
-Line(21) = {1, 11};
-Line(22) = {2, 12};
-Line(23) = {3, 13};
-Line(24) = {4, 14};
-
-// tip
-Ellipse(31) = {11, 10, 10, 21};
-Ellipse(32) = {12, 10, 10, 21};
-Ellipse(33) = {13, 10, 10, 21};
-Ellipse(34) = {14, 10, 10, 21};
-
-// box
-Line(41) = {41, 31};
-Line(42) = {31, 32};
-Line(43) = {32, 49};
-Line(44) = {49, 33};
-Line(45) = {33, 34};
-Line(46) = {34, 41};
-Line(47) = {43, 35};
-Line(48) = {35, 36};
-Line(49) = {36, 47};
-Line(50) = {47, 37};
-Line(51) = {37, 38};
-Line(52) = {38, 43};
-Line(53) = {31, 35};
-Line(54) = {32, 36};
-Line(55) = {33, 37};
-Line(56) = {34, 38};
-
-// midplane
-Line(61) = {1, 41};
-Line(62) = {11, 42};
-Line(63) = {11, 44};
-Line(64) = {21, 45};
-Line(65) = {13, 46};
-Line(66) = {13, 48};
-Line(67) = {3, 49};
-Line(68) = {41, 42};
-Line(69) = {42, 43};
-Line(70) = {43, 44};
-Line(71) = {44, 45};
-Line(72) = {45, 46};
-Line(73) = {46, 47};
-Line(74) = {47, 48};
-Line(75) = {48, 49};
-
-// cylinder
-Line Loop(1) = {1, 22, -11, -21};
-Surface(1) = {-1};
-Line Loop(2) = {2, 23, -12, -22};
-Surface(2) = {-2};
-Line Loop(3) = {3, 24, -13, -23};
-Surface(3) = {-3};
-Line Loop(4) = {4, 21, -14, -24};
-Surface(4) = {-4};
-//tip
-Line Loop(5) = {11, 32, -31};
-Surface(5) = {-5};
-Line Loop(6) = {12, 33, -32};
-Surface(6) = {-6};
-Line Loop(7) = {13, 34, -33};
-Surface(7) = {-7};
-Line Loop(8) = {14, 31, -34};
-Surface(8) = {-8};
-
-// symmetry
-Line Loop(9) = {61, 41, 42, 43, -67, -2, -1};
-Plane Surface(9) = {-9};
-Line Loop(10) = {4, 61, -46, -45, -44, -67, 3};
-Plane Surface(10) = {10};
-
-// upstream
-Line Loop(11) = {43, -75, -74, -49, -54};
-Plane Surface(11) = {11};
-Line Loop(12) = {75, 44, 55, -50, 74};
-Plane Surface(12) = {12};
-
-// downstream
-Line Loop(13) = {41, 53, -47, -69, -68};
-Plane Surface(13) = {13};
-Line Loop(14) = {68, 69, -52, -56, 46};
-Plane Surface(14) = {14};
-
-// farfield
-Line Loop(15) = {42, 54, -48, -53};
-Plane Surface(15) = {15};
-Line Loop(16) = {47, 48, 49, -73, -72, -71, -70};
-Plane Surface(16) = {16};
-Line Loop(17) = {51, 52, 70, 71, 72, 73, 50};
-Plane Surface(17) = {17};
-Line Loop(18) = {45, 56, -51, -55};
-Plane Surface(18) = {18};
-
-// wake
-Line Loop(19) = {61, 68, -62, -21};
-Surface(19) = {19};
-
-// internal
-Line Loop(20) = {62, 69, 70, -63};
-Surface(20) = {20};
-Line Loop(21) = {64, -71, -63, 31};
-Surface(21) = {-21};
-Line Loop(22) = {65, -72, -64, -33};
-Surface(22) = {-22};
-Line Loop(23) = {74, -66, 65, 73};
-Surface(23) = {23};
-Line Loop(24) = {75, -67, 23, 66};
-Surface(24) = {24};
-
-// field
-Surface Loop(1) = {15, 9, 13, 16, 11, 2, 6, 5, 1, 24, 23, 22, 21, 20, 19};
-Volume(1) = {1};
-Surface Loop(2) = {10, 4, 8, 7, 3, 14, 17, 18, 12, 24, 23, 22, 21, 20, 19};
-Volume(2) = {2};
-
-// Mesh
-// cylinder, farfield, symmetry, midplane
-MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24} = msh2D;
-// 3D
-Mesh.Algorithm3D = msh3D;
-Mesh.OptimizeNetgen = 1;
-Mesh.Smoothing = 10;
-Mesh.SmoothNormals = 1;
-
-// Physical
-Physical Line("teTip") = {21, 62};
-Physical Line("wakeTip") = {62};
-Physical Surface("cylinder") = {1, 2};
-Physical Surface("cylinder_") = {3, 4};
-Physical Surface("tip") = {5, 6, 7, 8};
-Physical Surface("symmetry") = {9};
-Physical Surface("symmetry_") = {10};
-Physical Surface("upstream") = {11, 12};
-Physical Surface("downstream") = {13};
-Physical Surface("downstream_") = {14};
-Physical Surface("farfield") = {15, 16, 17, 18};
-Physical Surface("wake") = {19};
-Physical Volume("field") = {1};
-Physical Volume("field_") = {2};
-
diff --git a/dart/models/n0012_3.geo b/dart/models/n0012_3.geo
deleted file mode 100644
index c8b02ca..0000000
--- a/dart/models/n0012_3.geo
+++ /dev/null
@@ -1,442 +0,0 @@
-/* Rectangular NACA0012 wing */
-// Initially generated by unsRgridWingGen.m
-// For gmsh 4, use line 334 instead of line 335 (Bezier <- BSpline)
-
-// Parameters
-// domain and mesh
-DefineConstant[ spn = { 2.0, Name "Wing span" } ];
-DefineConstant[ lgt = { 3.0, Name "Channel length" } ];
-DefineConstant[ wdt = { 3.0, Name "Channel width" } ];
-DefineConstant[ hgt = { 3.0, Name "Channel height" } ];
-DefineConstant[ msLe = { 0.05, Name "Leading edge mesh size" } ];
-DefineConstant[ msTe = { 0.05, Name "Trailing edge mesh size" } ];
-DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ];
-
-//// GEOMETRY
-
-
-/// Points
-
-// Airfoil 1: naca0012, 101 points
-Point(1) = {1.000000,0.000000,0.000000,msTe};
-Point(2) = {0.999013,0.000000,0.000143};
-Point(3) = {0.996057,0.000000,0.000572};
-Point(4) = {0.991144,0.000000,0.001280};
-Point(5) = {0.984292,0.000000,0.002260};
-Point(6) = {0.975528,0.000000,0.003501};
-Point(7) = {0.964888,0.000000,0.004990};
-Point(8) = {0.952414,0.000000,0.006710};
-Point(9) = {0.938153,0.000000,0.008643};
-Point(10) = {0.922164,0.000000,0.010770};
-Point(11) = {0.904508,0.000000,0.013071,msTe};
-Point(12) = {0.885257,0.000000,0.015523};
-Point(13) = {0.864484,0.000000,0.018106};
-Point(14) = {0.842274,0.000000,0.020795};
-Point(15) = {0.818712,0.000000,0.023569};
-Point(16) = {0.793893,0.000000,0.026405};
-Point(17) = {0.767913,0.000000,0.029279};
-Point(18) = {0.740877,0.000000,0.032168};
-Point(19) = {0.712890,0.000000,0.035048};
-Point(20) = {0.684062,0.000000,0.037896};
-Point(21) = {0.654508,0.000000,0.040686};
-Point(22) = {0.624345,0.000000,0.043394};
-Point(23) = {0.593691,0.000000,0.045992};
-Point(24) = {0.562667,0.000000,0.048455};
-Point(25) = {0.531395,0.000000,0.050754};
-Point(26) = {0.500000,0.000000,0.052862};
-Point(27) = {0.468605,0.000000,0.054749};
-Point(28) = {0.437333,0.000000,0.056390};
-Point(29) = {0.406309,0.000000,0.057755};
-Point(30) = {0.375655,0.000000,0.058819};
-Point(31) = {0.345492,0.000000,0.059557};
-Point(32) = {0.315938,0.000000,0.059947};
-Point(33) = {0.287110,0.000000,0.059971,msTe};
-Point(34) = {0.259123,0.000000,0.059614};
-Point(35) = {0.232087,0.000000,0.058863};
-Point(36) = {0.206107,0.000000,0.057712};
-Point(37) = {0.181288,0.000000,0.056159};
-Point(38) = {0.157726,0.000000,0.054206};
-Point(39) = {0.135516,0.000000,0.051862};
-Point(40) = {0.114743,0.000000,0.049138};
-Point(41) = {0.095492,0.000000,0.046049};
-Point(42) = {0.077836,0.000000,0.042615};
-Point(43) = {0.061847,0.000000,0.038859};
-Point(44) = {0.047586,0.000000,0.034803};
-Point(45) = {0.035112,0.000000,0.030473};
-Point(46) = {0.024472,0.000000,0.025893};
-Point(47) = {0.015708,0.000000,0.021088};
-Point(48) = {0.008856,0.000000,0.016078};
-Point(49) = {0.003943,0.000000,0.010884};
-Point(50) = {0.000987,0.000000,0.005521};
-Point(51) = {0.000000,0.000000,0.000000,msLe};
-Point(52) = {0.000987,0.000000,-0.005521};
-Point(53) = {0.003943,0.000000,-0.010884};
-Point(54) = {0.008856,0.000000,-0.016078};
-Point(55) = {0.015708,0.000000,-0.021088};
-Point(56) = {0.024472,0.000000,-0.025893};
-Point(57) = {0.035112,0.000000,-0.030473};
-Point(58) = {0.047586,0.000000,-0.034803};
-Point(59) = {0.061847,0.000000,-0.038859};
-Point(60) = {0.077836,0.000000,-0.042615};
-Point(61) = {0.095492,0.000000,-0.046049};
-Point(62) = {0.114743,0.000000,-0.049138};
-Point(63) = {0.135516,0.000000,-0.051862};
-Point(64) = {0.157726,0.000000,-0.054206};
-Point(65) = {0.181288,0.000000,-0.056159};
-Point(66) = {0.206107,0.000000,-0.057712};
-Point(67) = {0.232087,0.000000,-0.058863};
-Point(68) = {0.259123,0.000000,-0.059614};
-Point(69) = {0.287110,0.000000,-0.059971,msTe};
-Point(70) = {0.315938,0.000000,-0.059947};
-Point(71) = {0.345492,0.000000,-0.059557};
-Point(72) = {0.375655,0.000000,-0.058819};
-Point(73) = {0.406309,0.000000,-0.057755};
-Point(74) = {0.437333,0.000000,-0.056390};
-Point(75) = {0.468605,0.000000,-0.054749};
-Point(76) = {0.500000,0.000000,-0.052862};
-Point(77) = {0.531395,0.000000,-0.050754};
-Point(78) = {0.562667,0.000000,-0.048455};
-Point(79) = {0.593691,0.000000,-0.045992};
-Point(80) = {0.624345,0.000000,-0.043394};
-Point(81) = {0.654508,0.000000,-0.040686};
-Point(82) = {0.684062,0.000000,-0.037896};
-Point(83) = {0.712890,0.000000,-0.035048};
-Point(84) = {0.740877,0.000000,-0.032168};
-Point(85) = {0.767913,0.000000,-0.029279};
-Point(86) = {0.793893,0.000000,-0.026405};
-Point(87) = {0.818712,0.000000,-0.023569};
-Point(88) = {0.842274,0.000000,-0.020795};
-Point(89) = {0.864484,0.000000,-0.018106};
-Point(90) = {0.885257,0.000000,-0.015523};
-Point(91) = {0.904508,0.000000,-0.013071,msTe};
-Point(92) = {0.922164,0.000000,-0.010770};
-Point(93) = {0.938153,0.000000,-0.008643};
-Point(94) = {0.952414,0.000000,-0.006710};
-Point(95) = {0.964888,0.000000,-0.004990};
-Point(96) = {0.975528,0.000000,-0.003501};
-Point(97) = {0.984292,0.000000,-0.002260};
-Point(98) = {0.991144,0.000000,-0.001280};
-Point(99) = {0.996057,0.000000,-0.000572};
-Point(100) = {0.999013,0.000000,-0.000143,msTe};
-
-// Airfoil 2: naca0012, 101 points
-Point(102) = {1.000000,spn,0.000000,msTe};
-Point(103) = {0.999013,spn,0.000143};
-Point(104) = {0.996057,spn,0.000572};
-Point(105) = {0.991144,spn,0.001280};
-Point(106) = {0.984292,spn,0.002260};
-Point(107) = {0.975528,spn,0.003501};
-Point(108) = {0.964888,spn,0.004990};
-Point(109) = {0.952414,spn,0.006710};
-Point(110) = {0.938153,spn,0.008643};
-Point(111) = {0.922164,spn,0.010770};
-Point(112) = {0.904508,spn,0.013071,msTe};
-Point(113) = {0.885257,spn,0.015523};
-Point(114) = {0.864484,spn,0.018106};
-Point(115) = {0.842274,spn,0.020795};
-Point(116) = {0.818712,spn,0.023569};
-Point(117) = {0.793893,spn,0.026405};
-Point(118) = {0.767913,spn,0.029279};
-Point(119) = {0.740877,spn,0.032168};
-Point(120) = {0.712890,spn,0.035048};
-Point(121) = {0.684062,spn,0.037896};
-Point(122) = {0.654508,spn,0.040686};
-Point(123) = {0.624345,spn,0.043394};
-Point(124) = {0.593691,spn,0.045992};
-Point(125) = {0.562667,spn,0.048455};
-Point(126) = {0.531395,spn,0.050754};
-Point(127) = {0.500000,spn,0.052862};
-Point(128) = {0.468605,spn,0.054749};
-Point(129) = {0.437333,spn,0.056390};
-Point(130) = {0.406309,spn,0.057755};
-Point(131) = {0.375655,spn,0.058819};
-Point(132) = {0.345492,spn,0.059557};
-Point(133) = {0.315938,spn,0.059947};
-Point(134) = {0.287110,spn,0.059971,msTe};
-Point(135) = {0.259123,spn,0.059614};
-Point(136) = {0.232087,spn,0.058863};
-Point(137) = {0.206107,spn,0.057712};
-Point(138) = {0.181288,spn,0.056159};
-Point(139) = {0.157726,spn,0.054206};
-Point(140) = {0.135516,spn,0.051862};
-Point(141) = {0.114743,spn,0.049138};
-Point(142) = {0.095492,spn,0.046049};
-Point(143) = {0.077836,spn,0.042615};
-Point(144) = {0.061847,spn,0.038859};
-Point(145) = {0.047586,spn,0.034803};
-Point(146) = {0.035112,spn,0.030473};
-Point(147) = {0.024472,spn,0.025893};
-Point(148) = {0.015708,spn,0.021088};
-Point(149) = {0.008856,spn,0.016078};
-Point(150) = {0.003943,spn,0.010884};
-Point(151) = {0.000987,spn,0.005521};
-Point(152) = {0.000000,spn,0.000000,msLe};
-Point(153) = {0.000987,spn,-0.005521};
-Point(154) = {0.003943,spn,-0.010884};
-Point(155) = {0.008856,spn,-0.016078};
-Point(156) = {0.015708,spn,-0.021088};
-Point(157) = {0.024472,spn,-0.025893};
-Point(158) = {0.035112,spn,-0.030473};
-Point(159) = {0.047586,spn,-0.034803};
-Point(160) = {0.061847,spn,-0.038859};
-Point(161) = {0.077836,spn,-0.042615};
-Point(162) = {0.095492,spn,-0.046049};
-Point(163) = {0.114743,spn,-0.049138};
-Point(164) = {0.135516,spn,-0.051862};
-Point(165) = {0.157726,spn,-0.054206};
-Point(166) = {0.181288,spn,-0.056159};
-Point(167) = {0.206107,spn,-0.057712};
-Point(168) = {0.232087,spn,-0.058863};
-Point(169) = {0.259123,spn,-0.059614};
-Point(170) = {0.287110,spn,-0.059971,msTe};
-Point(171) = {0.315938,spn,-0.059947};
-Point(172) = {0.345492,spn,-0.059557};
-Point(173) = {0.375655,spn,-0.058819};
-Point(174) = {0.406309,spn,-0.057755};
-Point(175) = {0.437333,spn,-0.056390};
-Point(176) = {0.468605,spn,-0.054749};
-Point(177) = {0.500000,spn,-0.052862};
-Point(178) = {0.531395,spn,-0.050754};
-Point(179) = {0.562667,spn,-0.048455};
-Point(180) = {0.593691,spn,-0.045992};
-Point(181) = {0.624345,spn,-0.043394};
-Point(182) = {0.654508,spn,-0.040686};
-Point(183) = {0.684062,spn,-0.037896};
-Point(184) = {0.712890,spn,-0.035048};
-Point(185) = {0.740877,spn,-0.032168};
-Point(186) = {0.767913,spn,-0.029279};
-Point(187) = {0.793893,spn,-0.026405};
-Point(188) = {0.818712,spn,-0.023569};
-Point(189) = {0.842274,spn,-0.020795};
-Point(190) = {0.864484,spn,-0.018106};
-Point(191) = {0.885257,spn,-0.015523};
-Point(192) = {0.904508,spn,-0.013071,msTe};
-Point(193) = {0.922164,spn,-0.010770};
-Point(194) = {0.938153,spn,-0.008643};
-Point(195) = {0.952414,spn,-0.006710};
-Point(196) = {0.964888,spn,-0.004990};
-Point(197) = {0.975528,spn,-0.003501};
-Point(198) = {0.984292,spn,-0.002260};
-Point(199) = {0.991144,spn,-0.001280};
-Point(200) = {0.996057,spn,-0.000572};
-Point(201) = {0.999013,spn,-0.000143,msTe};
-
-// Tip:
-Point(5101) = {0.999013,spn,0.000000};
-Point(5102) = {0.996057,spn+0.000125,0.000000};
-Point(5103) = {0.991144,spn+0.000331,0.000000};
-Point(5104) = {0.984292,spn+0.000620,0.000000};
-Point(5105) = {0.975528,spn+0.000989,0.000000};
-Point(5106) = {0.964888,spn+0.001437,0.000000};
-Point(5107) = {0.952414,spn+0.001963,0.000000};
-Point(5108) = {0.938153,spn+0.002563,0.000000};
-Point(5109) = {0.922164,spn+0.003237,0.000000};
-Point(5110) = {0.904508,spn+0.003981,0.000000,msTe};
-Point(5111) = {0.885257,spn+0.004791,0.000000};
-Point(5112) = {0.864484,spn+0.005666,0.000000};
-Point(5113) = {0.842274,spn+0.006602,0.000000};
-Point(5114) = {0.818712,spn+0.007594,0.000000};
-Point(5115) = {0.793893,spn+0.008640,0.000000};
-Point(5116) = {0.767913,spn+0.009734,0.000000};
-Point(5117) = {0.740877,spn+0.010873,0.000000};
-Point(5118) = {0.712890,spn+0.012052,0.000000};
-Point(5119) = {0.684062,spn+0.013266,0.000000};
-Point(5120) = {0.654508,spn+0.014511,0.000000};
-Point(5121) = {0.624345,spn+0.015781,0.000000};
-Point(5122) = {0.593691,spn+0.017072,0.000000};
-Point(5123) = {0.562667,spn+0.018379,0.000000};
-Point(5124) = {0.531395,spn+0.019696,0.000000};
-Point(5125) = {0.500000,spn+0.021019,0.000000};
-Point(5126) = {0.468605,spn+0.022341,0.000000};
-Point(5127) = {0.437333,spn+0.023658,0.000000};
-Point(5128) = {0.406309,spn+0.024965,0.000000};
-Point(5129) = {0.375655,spn+0.026256,0.000000};
-Point(5130) = {0.345492,spn+0.027526,0.000000};
-Point(5131) = {0.315938,spn+0.028771,0.000000};
-Point(5132) = {0.287110,spn+0.029985,0.000000,msTe};
-Point(5133) = {0.019492,spn+0.041801,0.000000};
-
-// Dummy tip center:
-Point(5349) = {0.904508,spn,0.000000};
-Point(5350) = {0.287110,spn,0.000000};
-
-// Box
-Point(10001) = {1+lgt, 0, 0, msF};
-Point(10002) = {1+lgt, 0, hgt, msF};
-Point(10003) = {-lgt, 0, hgt, msF};
-Point(10004) = {-lgt, 0, -hgt, msF};
-Point(10005) = {1+lgt, 0, -hgt, msF};
-Point(10011) = {1+lgt, wdt, hgt, msF};
-Point(10012) = {-lgt, wdt, hgt, msF};
-Point(10013) = {-lgt, wdt, -hgt, msF};
-Point(10014) = {1+lgt, wdt, -hgt, msF};
-
-// Wake
-Point(10021) = {1+lgt, spn, 0, msF};
-
-/// Lines
-
-// Airfoil 1:
-Spline(1) = {1,2,3,4,5,6,7,8,9,10,11};
-Spline(2) = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33};
-Spline(3) = {33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51};
-Spline(4) = {51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69};
-Spline(5) = {69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91};
-Spline(6) = {91,92,93,94,95,96,97,98,99,100,1};
-
-// Airfoil 2:
-Spline(7) = {102,103,104,105,106,107,108,109,110,111,112};
-Spline(8) = {112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134};
-Spline(9) = {134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152};
-Spline(10) = {152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170};
-Spline(11) = {170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192};
-Spline(12) = {192,193,194,195,196,197,198,199,200,201,102};
-
-// Airfoil 1 to airfoil 2:
-Line(61) = {1,102};
-Line(62) = {11,112};
-Line(63) = {33,134};
-Line(64) = {51,152};
-Line(65) = {69,170};
-Line(66) = {91,192};
-
-// Tip:
-Spline(121) = {102,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110};
-Spline(122) = {5110,5111,5112,5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132};
-If(GMSH_MAJOR_VERSION >= 4)
- Bezier(123) = {5132,5133,152};
-Else
- BSpline(123) = {5132,5133,152};
-EndIf
-Ellipse(124) = {112,5349,112,5110};
-Ellipse(125) = {134,5350,134,5132};
-Ellipse(126) = {170,5350,170,5132};
-Ellipse(127) = {192,5349,192,5110};
-
-// Box
-Line(201) = {10001,10002};
-Line(202) = {10002,10003};
-Line(203) = {10003,10004};
-Line(204) = {10004,10005};
-Line(205) = {10005,10001};
-Line(211) = {10011,10012};
-Line(212) = {10012,10013};
-Line(213) = {10013,10014};
-Line(214) = {10014,10011};
-Line(221) = {1, 10001};
-Line(222) = {102, 10021};
-Line(231) = {10002,10011};
-Line(232) = {10003,10012};
-Line(233) = {10004,10013};
-Line(234) = {10005,10014};
-
-// Wake
-Line(241) = {10001, 10021};
-
-/// Line loops & Surfaces
-
-// Wing 1:
-Line Loop(1) = {1,62,-7,-61};
-Line Loop(2) = {2,63,-8,-62};
-Line Loop(3) = {3,64,-9,-63};
-Line Loop(4) = {4,65,-10,-64};
-Line Loop(5) = {5,66,-11,-65};
-Line Loop(6) = {6,61,-12,-66};
-Surface(1) = {-1};
-Surface(2) = {-2};
-Surface(3) = {-3};
-Surface(4) = {-4};
-Surface(5) = {-5};
-Surface(6) = {-6};
-
-// Wingtip:
-Line Loop(11) = {7,124,-121};
-Line Loop(12) = {8,125,-122,-124};
-Line Loop(13) = {9,-123,-125};
-Line Loop(14) = {10,126,123};
-Line Loop(15) = {11,127,122,-126};
-Line Loop(16) = {12,121,-127};
-Surface(11) = {-11};
-Surface(12) = {-12};
-Surface(13) = {-13};
-Surface(14) = {-14};
-Surface(15) = {-15};
-Surface(16) = {-16};
-
-// Symmetry
-Line Loop(21) = {201, 202, 203, 204, 205};
-Line Loop(22) = {1, 2, 3, 4, 5, 6};
-Plane Surface(21) = {-21, 22};
-
-// Downsteam
-Line Loop(31) = {201, 231, -214, -234, 205};
-Plane Surface(31) = {31};
-
-// Farfield
-Line Loop(41) = {233, -212, -232, 203};
-Plane Surface(41) = {41};
-Line Loop(42) = {204, 234, -213, -233};
-Plane Surface(42) = {42};
-Line Loop(43) = {202, 232, -211, -231};
-Plane Surface(43) = {43};
-Line Loop(44) = {213, 214, 211, 212};
-Plane Surface(44) = {44};
-
-// Wake
-Line Loop(51) = {221, 241, -222, -61};
-Surface(51) = {51};
-
-/// Volume
-Surface Loop(1) = {1,2,3,4,5,6,11,12,13,14,15,16,21,31,41,42,43,44};
-Volume(1) = {1};
-
-/// Embbeded
-Line{221} In Surface{21};
-Line{241} In Surface{31};
-Surface{51} In Volume{1};
-
-//// MESHING ALGORITHM
-
-/// 2D:
-
-///Wing, farfield and symmetry plane:
-MeshAlgorithm Surface{1,2,3,4,5,6,21,31,41,42,43,44,51} = 5;
-
-///Tip:
-MeshAlgorithm Surface{11,12,13,14,15,16} = 1;
-
-/// 3D:
-
-Mesh.Algorithm3D = 2;
-//Mesh.OptimizeNetgen = 1;
-Mesh.Optimize = 1;
-Mesh.Smoothing = 10;
-Mesh.SmoothNormals = 1;
-
-
-
-//// PHYSICAL GROUPS
-
-// Trailing edge and wake tip
-Physical Line("wakeTip") = {222};
-Physical Line("teTip") = {61, 222};
-
-// Internal Field:
-Physical Volume("field") = {1};
-
-// Wing:
-Physical Surface("wing") = {1,2,3,11,12,13};
-Physical Surface("wing_") = {4,5,6,14,15,16};
-
-// Symmetry:
-Physical Surface("symmetry") = {21};
-
-// Downstream:
-Physical Surface("downstream") = {31};
-
-// Farfield:
-Physical Surface("upstream") = {41};
-Physical Surface("farfield") = {42,43,44};
-
-// Wake:
-Physical Surface("wake") = {51};
diff --git a/dart/models/oneraM6.geo b/dart/models/oneraM6.geo
index 6c14627..25b777b 100644
--- a/dart/models/oneraM6.geo
+++ b/dart/models/oneraM6.geo
@@ -559,7 +559,7 @@ Mesh.SmoothNormals = 1;
// Trailing edge and wake tip
Physical Line("wakeTip") = {162};
-Physical Line("teTip") = {61, 162};
+Physical Line("te") = {61};
// Internal Field:
Physical Volume("field") = {1};
diff --git a/dart/models/rae_25.geo b/dart/models/rae_25.geo
new file mode 100644
index 0000000..1deee3c
--- /dev/null
+++ b/dart/models/rae_25.geo
@@ -0,0 +1,400 @@
+/* 2.5D extruded RAE 2822 wing */
+
+// Parameters
+// domain
+DefineConstant[ spn = { 1.0, Name "Wing span" },
+ lgt = { 6.0, Name "Channel length" },
+ hgt = { 6.0, Name "Channel height" },
+ msN = { 0.05, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" },
+ msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ];
+
+// Geometry
+// wing
+Point(1) = {1.000000, -spn/2, 0.000000, msN};
+Point(2) = {0.999398, -spn/2, 0.000128, msN};
+Point(3) = {0.997592, -spn/2, 0.000510, msN};
+Point(4) = {0.994588, -spn/2, 0.001137, msN};
+Point(5) = {0.990393, -spn/2, 0.002001, msN};
+Point(6) = {0.985016, -spn/2, 0.003092, msN};
+Point(7) = {0.978470, -spn/2, 0.004401, msN};
+Point(8) = {0.970772, -spn/2, 0.005915, msN};
+Point(9) = {0.961940, -spn/2, 0.007622, msN};
+Point(10) = {0.951995, -spn/2, 0.009508, msN};
+Point(11) = {0.940961, -spn/2, 0.011562, msN};
+Point(12) = {0.928864, -spn/2, 0.013769, msN};
+Point(13) = {0.915735, -spn/2, 0.016113, msN};
+Point(14) = {0.901604, -spn/2, 0.018580, msN};
+Point(15) = {0.886505, -spn/2, 0.021153, msN};
+Point(16) = {0.870476, -spn/2, 0.023817, msN};
+Point(17) = {0.853553, -spn/2, 0.026554, msN};
+Point(18) = {0.835779, -spn/2, 0.029347, msN};
+Point(19) = {0.817197, -spn/2, 0.032176, msN};
+Point(20) = {0.797850, -spn/2, 0.035017, msN};
+Point(21) = {0.777785, -spn/2, 0.037847, msN};
+Point(22) = {0.757051, -spn/2, 0.040641, msN};
+Point(23) = {0.735698, -spn/2, 0.043377, msN};
+Point(24) = {0.713778, -spn/2, 0.046029, msN};
+Point(25) = {0.691342, -spn/2, 0.048575, msN};
+Point(26) = {0.668445, -spn/2, 0.050993, msN};
+Point(27) = {0.645142, -spn/2, 0.053258, msN};
+Point(28) = {0.621490, -spn/2, 0.055344, msN};
+Point(29) = {0.597545, -spn/2, 0.057218, msN};
+Point(30) = {0.573365, -spn/2, 0.058845, msN};
+Point(31) = {0.549009, -spn/2, 0.060194, msN};
+Point(32) = {0.524534, -spn/2, 0.061254, msN};
+Point(33) = {0.500000, -spn/2, 0.062029, msN};
+Point(34) = {0.475466, -spn/2, 0.062530, msN};
+Point(35) = {0.450991, -spn/2, 0.062774, msN};
+Point(36) = {0.426635, -spn/2, 0.062779, msN};
+Point(37) = {0.402455, -spn/2, 0.062562, msN};
+Point(38) = {0.378510, -spn/2, 0.062133, msN};
+Point(39) = {0.354858, -spn/2, 0.061497, msN};
+Point(40) = {0.331555, -spn/2, 0.060660, msN};
+Point(41) = {0.308658, -spn/2, 0.059629, msN};
+Point(42) = {0.286222, -spn/2, 0.058414, msN};
+Point(43) = {0.264302, -spn/2, 0.057026, msN};
+Point(44) = {0.242949, -spn/2, 0.055470, msN};
+Point(45) = {0.222215, -spn/2, 0.053753, msN};
+Point(46) = {0.202150, -spn/2, 0.051885, msN};
+Point(47) = {0.182803, -spn/2, 0.049874, msN};
+Point(48) = {0.164221, -spn/2, 0.047729, msN};
+Point(49) = {0.146447, -spn/2, 0.045457, msN};
+Point(50) = {0.129524, -spn/2, 0.043071, msN};
+Point(51) = {0.113495, -spn/2, 0.040585, msN};
+Point(52) = {0.098396, -spn/2, 0.038011, msN};
+Point(53) = {0.084265, -spn/2, 0.035360, msN};
+Point(54) = {0.071136, -spn/2, 0.032644, msN};
+Point(55) = {0.059039, -spn/2, 0.029874, msN};
+Point(56) = {0.048005, -spn/2, 0.027062, msN};
+Point(57) = {0.038060, -spn/2, 0.024219, msN};
+Point(58) = {0.029228, -spn/2, 0.021348, msN};
+Point(59) = {0.021530, -spn/2, 0.018441, msN};
+Point(60) = {0.014984, -spn/2, 0.015489, msN};
+Point(61) = {0.009607, -spn/2, 0.012480, msN};
+Point(62) = {0.005412, -spn/2, 0.009416, msN};
+Point(63) = {0.002408, -spn/2, 0.006306, msN};
+Point(64) = {0.000602, -spn/2, 0.003165, msN};
+Point(65) = {0.000000, -spn/2, 0.000000, msN};
+Point(66) = {0.000602, -spn/2, -0.003160, msN};
+Point(67) = {0.002408, -spn/2, -0.006308, msN};
+Point(68) = {0.005412, -spn/2, -0.009443, msN};
+Point(69) = {0.009607, -spn/2, -0.012559, msN};
+Point(70) = {0.014984, -spn/2, -0.015649, msN};
+Point(71) = {0.021530, -spn/2, -0.018707, msN};
+Point(72) = {0.029228, -spn/2, -0.021722, msN};
+Point(73) = {0.038060, -spn/2, -0.024685, msN};
+Point(74) = {0.048005, -spn/2, -0.027586, msN};
+Point(75) = {0.059039, -spn/2, -0.030416, msN};
+Point(76) = {0.071136, -spn/2, -0.033170, msN};
+Point(77) = {0.084265, -spn/2, -0.035843, msN};
+Point(78) = {0.098396, -spn/2, -0.038431, msN};
+Point(79) = {0.113495, -spn/2, -0.040929, msN};
+Point(80) = {0.129524, -spn/2, -0.043326, msN};
+Point(81) = {0.146447, -spn/2, -0.045610, msN};
+Point(82) = {0.164221, -spn/2, -0.047773, msN};
+Point(83) = {0.182803, -spn/2, -0.049805, msN};
+Point(84) = {0.202150, -spn/2, -0.051694, msN};
+Point(85) = {0.222215, -spn/2, -0.053427, msN};
+Point(86) = {0.242949, -spn/2, -0.054994, msN};
+Point(87) = {0.264302, -spn/2, -0.056376, msN};
+Point(88) = {0.286222, -spn/2, -0.057547, msN};
+Point(89) = {0.308658, -spn/2, -0.058459, msN};
+Point(90) = {0.331555, -spn/2, -0.059046, msN};
+Point(91) = {0.354858, -spn/2, -0.059236, msN};
+Point(92) = {0.378510, -spn/2, -0.058974, msN};
+Point(93) = {0.402455, -spn/2, -0.058224, msN};
+Point(94) = {0.426635, -spn/2, -0.056979, msN};
+Point(95) = {0.450991, -spn/2, -0.055257, msN};
+Point(96) = {0.475466, -spn/2, -0.053099, msN};
+Point(97) = {0.500000, -spn/2, -0.050563, msN};
+Point(98) = {0.524534, -spn/2, -0.047719, msN};
+Point(99) = {0.549009, -spn/2, -0.044642, msN};
+Point(100) = {0.573365, -spn/2, -0.041397, msN};
+Point(101) = {0.597545, -spn/2, -0.038043, msN};
+Point(102) = {0.621490, -spn/2, -0.034631, msN};
+Point(103) = {0.645142, -spn/2, -0.031207, msN};
+Point(104) = {0.668445, -spn/2, -0.027814, msN};
+Point(105) = {0.691342, -spn/2, -0.024495, msN};
+Point(106) = {0.713778, -spn/2, -0.021289, msN};
+Point(107) = {0.735698, -spn/2, -0.018232, msN};
+Point(108) = {0.757051, -spn/2, -0.015357, msN};
+Point(109) = {0.777785, -spn/2, -0.012690, msN};
+Point(110) = {0.797850, -spn/2, -0.010244, msN};
+Point(111) = {0.817197, -spn/2, -0.008027, msN};
+Point(112) = {0.835779, -spn/2, -0.006048, msN};
+Point(113) = {0.853553, -spn/2, -0.004314, msN};
+Point(114) = {0.870476, -spn/2, -0.002829, msN};
+Point(115) = {0.886505, -spn/2, -0.001592, msN};
+Point(116) = {0.901604, -spn/2, -0.000600, msN};
+Point(117) = {0.915735, -spn/2, 0.000157, msN};
+Point(118) = {0.928864, -spn/2, 0.000694, msN};
+Point(119) = {0.940961, -spn/2, 0.001033, msN};
+Point(120) = {0.951995, -spn/2, 0.001197, msN};
+Point(121) = {0.961940, -spn/2, 0.001212, msN};
+Point(122) = {0.970772, -spn/2, 0.001112, msN};
+Point(123) = {0.978470, -spn/2, 0.000935, msN};
+Point(124) = {0.985016, -spn/2, 0.000719, msN};
+Point(125) = {0.990393, -spn/2, 0.000497, msN};
+Point(126) = {0.994588, -spn/2, 0.000296, msN};
+Point(127) = {0.997592, -spn/2, 0.000137, msN};
+Point(128) = {0.999398, -spn/2, 0.000035, msN};
+
+Point(301) = {1.000000, spn/2, 0.000000, msN};
+Point(302) = {0.999398, spn/2, 0.000128, msN};
+Point(303) = {0.997592, spn/2, 0.000510, msN};
+Point(304) = {0.994588, spn/2, 0.001137, msN};
+Point(305) = {0.990393, spn/2, 0.002001, msN};
+Point(306) = {0.985016, spn/2, 0.003092, msN};
+Point(307) = {0.978470, spn/2, 0.004401, msN};
+Point(308) = {0.970772, spn/2, 0.005915, msN};
+Point(309) = {0.961940, spn/2, 0.007622, msN};
+Point(310) = {0.951995, spn/2, 0.009508, msN};
+Point(311) = {0.940961, spn/2, 0.011562, msN};
+Point(312) = {0.928864, spn/2, 0.013769, msN};
+Point(313) = {0.915735, spn/2, 0.016113, msN};
+Point(314) = {0.901604, spn/2, 0.018580, msN};
+Point(315) = {0.886505, spn/2, 0.021153, msN};
+Point(316) = {0.870476, spn/2, 0.023817, msN};
+Point(317) = {0.853553, spn/2, 0.026554, msN};
+Point(318) = {0.835779, spn/2, 0.029347, msN};
+Point(319) = {0.817197, spn/2, 0.032176, msN};
+Point(320) = {0.797850, spn/2, 0.035017, msN};
+Point(321) = {0.777785, spn/2, 0.037847, msN};
+Point(322) = {0.757051, spn/2, 0.040641, msN};
+Point(323) = {0.735698, spn/2, 0.043377, msN};
+Point(324) = {0.713778, spn/2, 0.046029, msN};
+Point(325) = {0.691342, spn/2, 0.048575, msN};
+Point(326) = {0.668445, spn/2, 0.050993, msN};
+Point(327) = {0.645142, spn/2, 0.053258, msN};
+Point(328) = {0.621490, spn/2, 0.055344, msN};
+Point(329) = {0.597545, spn/2, 0.057218, msN};
+Point(330) = {0.573365, spn/2, 0.058845, msN};
+Point(331) = {0.549009, spn/2, 0.060194, msN};
+Point(332) = {0.524534, spn/2, 0.061254, msN};
+Point(333) = {0.500000, spn/2, 0.062029, msN};
+Point(334) = {0.475466, spn/2, 0.062530, msN};
+Point(335) = {0.450991, spn/2, 0.062774, msN};
+Point(336) = {0.426635, spn/2, 0.062779, msN};
+Point(337) = {0.402455, spn/2, 0.062562, msN};
+Point(338) = {0.378510, spn/2, 0.062133, msN};
+Point(339) = {0.354858, spn/2, 0.061497, msN};
+Point(340) = {0.331555, spn/2, 0.060660, msN};
+Point(341) = {0.308658, spn/2, 0.059629, msN};
+Point(342) = {0.286222, spn/2, 0.058414, msN};
+Point(343) = {0.264302, spn/2, 0.057026, msN};
+Point(344) = {0.242949, spn/2, 0.055470, msN};
+Point(345) = {0.222215, spn/2, 0.053753, msN};
+Point(346) = {0.202150, spn/2, 0.051885, msN};
+Point(347) = {0.182803, spn/2, 0.049874, msN};
+Point(348) = {0.164221, spn/2, 0.047729, msN};
+Point(349) = {0.146447, spn/2, 0.045457, msN};
+Point(350) = {0.129524, spn/2, 0.043071, msN};
+Point(351) = {0.113495, spn/2, 0.040585, msN};
+Point(352) = {0.098396, spn/2, 0.038011, msN};
+Point(353) = {0.084265, spn/2, 0.035360, msN};
+Point(354) = {0.071136, spn/2, 0.032644, msN};
+Point(355) = {0.059039, spn/2, 0.029874, msN};
+Point(356) = {0.048005, spn/2, 0.027062, msN};
+Point(357) = {0.038060, spn/2, 0.024219, msN};
+Point(358) = {0.029228, spn/2, 0.021348, msN};
+Point(359) = {0.021530, spn/2, 0.018441, msN};
+Point(360) = {0.014984, spn/2, 0.015489, msN};
+Point(361) = {0.009607, spn/2, 0.012480, msN};
+Point(362) = {0.005412, spn/2, 0.009416, msN};
+Point(363) = {0.002408, spn/2, 0.006306, msN};
+Point(364) = {0.000602, spn/2, 0.003165, msN};
+Point(365) = {0.000000, spn/2, 0.000000, msN};
+Point(366) = {0.000602, spn/2, -0.003160, msN};
+Point(367) = {0.002408, spn/2, -0.006308, msN};
+Point(368) = {0.005412, spn/2, -0.009443, msN};
+Point(369) = {0.009607, spn/2, -0.012559, msN};
+Point(370) = {0.014984, spn/2, -0.015649, msN};
+Point(371) = {0.021530, spn/2, -0.018707, msN};
+Point(372) = {0.029228, spn/2, -0.021722, msN};
+Point(373) = {0.038060, spn/2, -0.024685, msN};
+Point(374) = {0.048005, spn/2, -0.027586, msN};
+Point(375) = {0.059039, spn/2, -0.030416, msN};
+Point(376) = {0.071136, spn/2, -0.033170, msN};
+Point(377) = {0.084265, spn/2, -0.035843, msN};
+Point(378) = {0.098396, spn/2, -0.038431, msN};
+Point(379) = {0.113495, spn/2, -0.040929, msN};
+Point(380) = {0.129524, spn/2, -0.043326, msN};
+Point(381) = {0.146447, spn/2, -0.045610, msN};
+Point(382) = {0.164221, spn/2, -0.047773, msN};
+Point(383) = {0.182803, spn/2, -0.049805, msN};
+Point(384) = {0.202150, spn/2, -0.051694, msN};
+Point(385) = {0.222215, spn/2, -0.053427, msN};
+Point(386) = {0.242949, spn/2, -0.054994, msN};
+Point(387) = {0.264302, spn/2, -0.056376, msN};
+Point(388) = {0.286222, spn/2, -0.057547, msN};
+Point(389) = {0.308658, spn/2, -0.058459, msN};
+Point(390) = {0.331555, spn/2, -0.059046, msN};
+Point(391) = {0.354858, spn/2, -0.059236, msN};
+Point(392) = {0.378510, spn/2, -0.058974, msN};
+Point(393) = {0.402455, spn/2, -0.058224, msN};
+Point(394) = {0.426635, spn/2, -0.056979, msN};
+Point(395) = {0.450991, spn/2, -0.055257, msN};
+Point(396) = {0.475466, spn/2, -0.053099, msN};
+Point(397) = {0.500000, spn/2, -0.050563, msN};
+Point(398) = {0.524534, spn/2, -0.047719, msN};
+Point(399) = {0.549009, spn/2, -0.044642, msN};
+Point(400) = {0.573365, spn/2, -0.041397, msN};
+Point(401) = {0.597545, spn/2, -0.038043, msN};
+Point(402) = {0.621490, spn/2, -0.034631, msN};
+Point(403) = {0.645142, spn/2, -0.031207, msN};
+Point(404) = {0.668445, spn/2, -0.027814, msN};
+Point(405) = {0.691342, spn/2, -0.024495, msN};
+Point(406) = {0.713778, spn/2, -0.021289, msN};
+Point(407) = {0.735698, spn/2, -0.018232, msN};
+Point(408) = {0.757051, spn/2, -0.015357, msN};
+Point(409) = {0.777785, spn/2, -0.012690, msN};
+Point(410) = {0.797850, spn/2, -0.010244, msN};
+Point(411) = {0.817197, spn/2, -0.008027, msN};
+Point(412) = {0.835779, spn/2, -0.006048, msN};
+Point(413) = {0.853553, spn/2, -0.004314, msN};
+Point(414) = {0.870476, spn/2, -0.002829, msN};
+Point(415) = {0.886505, spn/2, -0.001592, msN};
+Point(416) = {0.901604, spn/2, -0.000600, msN};
+Point(417) = {0.915735, spn/2, 0.000157, msN};
+Point(418) = {0.928864, spn/2, 0.000694, msN};
+Point(419) = {0.940961, spn/2, 0.001033, msN};
+Point(420) = {0.951995, spn/2, 0.001197, msN};
+Point(421) = {0.961940, spn/2, 0.001212, msN};
+Point(422) = {0.970772, spn/2, 0.001112, msN};
+Point(423) = {0.978470, spn/2, 0.000935, msN};
+Point(424) = {0.985016, spn/2, 0.000719, msN};
+Point(425) = {0.990393, spn/2, 0.000497, msN};
+Point(426) = {0.994588, spn/2, 0.000296, msN};
+Point(427) = {0.997592, spn/2, 0.000137, msN};
+Point(428) = {0.999398, spn/2, 0.000035, msN};
+
+// box
+Point(631) = {1+lgt/2, -spn/2, hgt/2, msF};
+Point(632) = {-lgt/2, -spn/2, hgt/2, msF};
+Point(633) = {-lgt/2, -spn/2, -hgt/2, msF};
+Point(634) = {1+lgt/2, -spn/2, -hgt/2, msF};
+Point(635) = {1+lgt/2, spn/2, hgt/2, msF};
+Point(636) = {-lgt/2, spn/2, hgt/2, msF};
+Point(637) = {-lgt/2, spn/2, -hgt/2, msF};
+Point(638) = {1+lgt/2, spn/2, -hgt/2, msF};
+
+// midplane
+Point(641) = {1+lgt/2, -spn/2, 0, msF};
+Point(642) = {1+lgt/2, spn/2, 0, msF};
+Point(643) = {-lgt/2, spn/2, 0, msF};
+Point(644) = {-lgt/2, -spn/2, 0, msF};
+
+// wing
+Spline(1) = {1:44};
+Spline(2) = {44:65};
+Spline(3) = {65:86};
+Spline(4) = {86:128,1};
+Spline(5) = {301:344};
+Spline(6) = {344:365};
+Spline(7) = {365:386};
+Spline(8) = {386:428,301};
+Line(9) = {1, 301};
+Line(10) = {44, 344};
+Line(11) = {65, 365};
+Line(12) = {86, 386};
+
+// box
+Line(21) = {641, 631};
+Line(22) = {631, 632};
+Line(23) = {632, 644};
+Line(24) = {644, 633};
+Line(25) = {633, 634};
+Line(26) = {634, 641};
+Line(27) = {642, 635};
+Line(28) = {635, 636};
+Line(29) = {636, 643};
+Line(30) = {643, 637};
+Line(31) = {637, 638};
+Line(32) = {638, 642};
+Line(33) = {631, 635};
+Line(34) = {632, 636};
+Line(35) = {633, 637};
+Line(36) = {634, 638};
+
+// midplane
+Line(41) = {1, 641};
+Line(42) = {301, 642};
+Line(43) = {365, 643};
+Line(44) = {65, 644};
+Line(45) = {641, 642};
+Line(46) = {643, 644};
+
+// wing
+Line Loop(1) = {1, 10, -5, -9};
+Surface(1) = {-1};
+Line Loop(2) = {2, 11, -6, -10};
+Surface(2) = {-2};
+Line Loop(3) = {3, 12, -7, -11};
+Surface(3) = {-3};
+Line Loop(4) = {4, 9, -8, -12};
+Surface(4) = {-4};
+
+// symmetry
+Line Loop(5) = {41, 21, 22, 23, -44, -2, -1};
+Plane Surface(5) = {-5};
+Line Loop(12) = {42, 27, 28, 29, -43, -5, -6};
+Plane Surface(12) = {12};
+Line Loop(6) = {44, 24, 25, 26, -41, -4, -3};
+Plane Surface(6) = {-6};
+Line Loop(13) = {43, 30, 31, 32, -42, -8, -7};
+Plane Surface(13) = {13};
+
+// upstream
+Line Loop(7) = {46, -23, 34, 29};
+Plane Surface(7) = {-7};
+Line Loop(8) = {46, 24, 35, -30};
+Plane Surface(8) = {8};
+
+// downstream
+Line Loop(9) = {45, 27, -33, -21};
+Plane Surface(9) = {-9};
+Line Loop(10) = {45, -32, -36, 26};
+Plane Surface(10) = {10};
+
+// farfield
+Line Loop(11) = {22, 34, -28, -33};
+Plane Surface(11) = {11};
+Line Loop(14) = {25, 36, -31, -35};
+Plane Surface(14) = {14};
+
+// wake
+Line Loop(15) = {41, 45, -42, -9};
+Plane Surface(15) = {15};
+
+// internal
+Line Loop(16) = {11, 43, 46, -44};
+Plane Surface(16) = {16};
+
+// field
+Surface Loop(1) = {7, 5, 9, 12, 11, 1, 2, 16, 15};
+Volume(1) = {1};
+Surface Loop(2) = {6, 8, 14, 10, 13, 3, 4, 16, 15};
+Volume(2) = {2};
+
+// Mesh
+// 2D
+MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} = 5;
+// 3D
+Mesh.Algorithm3D = 2;
+Mesh.OptimizeNetgen = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+// Physical
+Physical Line("te") = {9};
+Physical Surface("wing") = {1, 2};
+Physical Surface("wing_") = {3, 4};
+Physical Surface("symmetry") = {5, 12};
+Physical Surface("symmetry_") = {6, 13};
+Physical Surface("upstream") = {7, 8};
+Physical Surface("downstream") = {9};
+Physical Surface("downstream_") = {10};
+Physical Surface("farfield") = {11, 14};
+Physical Surface("wake") = {15};
+Physical Volume("field") = {1};
+Physical Volume("field_") = {2};
diff --git a/dart/models/rae_3.geo b/dart/models/rae_3.geo
new file mode 100644
index 0000000..2c76fac
--- /dev/null
+++ b/dart/models/rae_3.geo
@@ -0,0 +1,430 @@
+/* RAE 2822 rectangular wing */
+
+// Parameters
+// domain and mesh
+DefineConstant[ spn = { 1.0, Name "Wing span" },
+ lgt = { 6.0, Name "Channel length" },
+ wdt = { 3.0, Name "Channel width" },
+ hgt = { 6.0, Name "Channel height" },
+ msN = { 0.05, Name "Leading edge mesh size" },
+ msF = { 1.0, Name "Farfield mesh size" } ];
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1
+Point(1) = {1.000000, 0.0000, 0.000000, msN};
+Point(2) = {0.999398, 0.0000, 0.000128, msN};
+Point(3) = {0.997592, 0.0000, 0.000510, msN};
+Point(4) = {0.994588, 0.0000, 0.001137, msN};
+Point(5) = {0.990393, 0.0000, 0.002001, msN};
+Point(6) = {0.985016, 0.0000, 0.003092, msN};
+Point(7) = {0.978470, 0.0000, 0.004401, msN};
+Point(8) = {0.970772, 0.0000, 0.005915, msN};
+Point(9) = {0.961940, 0.0000, 0.007622, msN};
+Point(10) = {0.951995, 0.0000, 0.009508, msN};
+Point(11) = {0.940961, 0.0000, 0.011562, msN};
+Point(12) = {0.928864, 0.0000, 0.013769, msN};
+Point(13) = {0.915735, 0.0000, 0.016113, msN};
+Point(14) = {0.901604, 0.0000, 0.018580, msN};
+Point(15) = {0.886505, 0.0000, 0.021153, msN};
+Point(16) = {0.870476, 0.0000, 0.023817, msN};
+Point(17) = {0.853553, 0.0000, 0.026554, msN};
+Point(18) = {0.835779, 0.0000, 0.029347, msN};
+Point(19) = {0.817197, 0.0000, 0.032176, msN};
+Point(20) = {0.797850, 0.0000, 0.035017, msN};
+Point(21) = {0.777785, 0.0000, 0.037847, msN};
+Point(22) = {0.757051, 0.0000, 0.040641, msN};
+Point(23) = {0.735698, 0.0000, 0.043377, msN};
+Point(24) = {0.713778, 0.0000, 0.046029, msN};
+Point(25) = {0.691342, 0.0000, 0.048575, msN};
+Point(26) = {0.668445, 0.0000, 0.050993, msN};
+Point(27) = {0.645142, 0.0000, 0.053258, msN};
+Point(28) = {0.621490, 0.0000, 0.055344, msN};
+Point(29) = {0.597545, 0.0000, 0.057218, msN};
+Point(30) = {0.573365, 0.0000, 0.058845, msN};
+Point(31) = {0.549009, 0.0000, 0.060194, msN};
+Point(32) = {0.524534, 0.0000, 0.061254, msN};
+Point(33) = {0.500000, 0.0000, 0.062029, msN};
+Point(34) = {0.475466, 0.0000, 0.062530, msN};
+Point(35) = {0.450991, 0.0000, 0.062774, msN};
+Point(36) = {0.426635, 0.0000, 0.062779, msN};
+Point(37) = {0.402455, 0.0000, 0.062562, msN};
+Point(38) = {0.378510, 0.0000, 0.062133, msN};
+Point(39) = {0.354858, 0.0000, 0.061497, msN};
+Point(40) = {0.331555, 0.0000, 0.060660, msN};
+Point(41) = {0.308658, 0.0000, 0.059629, msN};
+Point(42) = {0.286222, 0.0000, 0.058414, msN};
+Point(43) = {0.264302, 0.0000, 0.057026, msN};
+Point(44) = {0.242949, 0.0000, 0.055470, msN};
+Point(45) = {0.222215, 0.0000, 0.053753, msN};
+Point(46) = {0.202150, 0.0000, 0.051885, msN};
+Point(47) = {0.182803, 0.0000, 0.049874, msN};
+Point(48) = {0.164221, 0.0000, 0.047729, msN};
+Point(49) = {0.146447, 0.0000, 0.045457, msN};
+Point(50) = {0.129524, 0.0000, 0.043071, msN};
+Point(51) = {0.113495, 0.0000, 0.040585, msN};
+Point(52) = {0.098396, 0.0000, 0.038011, msN};
+Point(53) = {0.084265, 0.0000, 0.035360, msN};
+Point(54) = {0.071136, 0.0000, 0.032644, msN};
+Point(55) = {0.059039, 0.0000, 0.029874, msN};
+Point(56) = {0.048005, 0.0000, 0.027062, msN};
+Point(57) = {0.038060, 0.0000, 0.024219, msN};
+Point(58) = {0.029228, 0.0000, 0.021348, msN};
+Point(59) = {0.021530, 0.0000, 0.018441, msN};
+Point(60) = {0.014984, 0.0000, 0.015489, msN};
+Point(61) = {0.009607, 0.0000, 0.012480, msN};
+Point(62) = {0.005412, 0.0000, 0.009416, msN};
+Point(63) = {0.002408, 0.0000, 0.006306, msN};
+Point(64) = {0.000602, 0.0000, 0.003165, msN};
+Point(65) = {0.000000, 0.0000, 0.000000, msN};
+Point(66) = {0.000602, 0.0000, -0.003160, msN};
+Point(67) = {0.002408, 0.0000, -0.006308, msN};
+Point(68) = {0.005412, 0.0000, -0.009443, msN};
+Point(69) = {0.009607, 0.0000, -0.012559, msN};
+Point(70) = {0.014984, 0.0000, -0.015649, msN};
+Point(71) = {0.021530, 0.0000, -0.018707, msN};
+Point(72) = {0.029228, 0.0000, -0.021722, msN};
+Point(73) = {0.038060, 0.0000, -0.024685, msN};
+Point(74) = {0.048005, 0.0000, -0.027586, msN};
+Point(75) = {0.059039, 0.0000, -0.030416, msN};
+Point(76) = {0.071136, 0.0000, -0.033170, msN};
+Point(77) = {0.084265, 0.0000, -0.035843, msN};
+Point(78) = {0.098396, 0.0000, -0.038431, msN};
+Point(79) = {0.113495, 0.0000, -0.040929, msN};
+Point(80) = {0.129524, 0.0000, -0.043326, msN};
+Point(81) = {0.146447, 0.0000, -0.045610, msN};
+Point(82) = {0.164221, 0.0000, -0.047773, msN};
+Point(83) = {0.182803, 0.0000, -0.049805, msN};
+Point(84) = {0.202150, 0.0000, -0.051694, msN};
+Point(85) = {0.222215, 0.0000, -0.053427, msN};
+Point(86) = {0.242949, 0.0000, -0.054994, msN};
+Point(87) = {0.264302, 0.0000, -0.056376, msN};
+Point(88) = {0.286222, 0.0000, -0.057547, msN};
+Point(89) = {0.308658, 0.0000, -0.058459, msN};
+Point(90) = {0.331555, 0.0000, -0.059046, msN};
+Point(91) = {0.354858, 0.0000, -0.059236, msN};
+Point(92) = {0.378510, 0.0000, -0.058974, msN};
+Point(93) = {0.402455, 0.0000, -0.058224, msN};
+Point(94) = {0.426635, 0.0000, -0.056979, msN};
+Point(95) = {0.450991, 0.0000, -0.055257, msN};
+Point(96) = {0.475466, 0.0000, -0.053099, msN};
+Point(97) = {0.500000, 0.0000, -0.050563, msN};
+Point(98) = {0.524534, 0.0000, -0.047719, msN};
+Point(99) = {0.549009, 0.0000, -0.044642, msN};
+Point(100) = {0.573365, 0.0000, -0.041397, msN};
+Point(101) = {0.597545, 0.0000, -0.038043, msN};
+Point(102) = {0.621490, 0.0000, -0.034631, msN};
+Point(103) = {0.645142, 0.0000, -0.031207, msN};
+Point(104) = {0.668445, 0.0000, -0.027814, msN};
+Point(105) = {0.691342, 0.0000, -0.024495, msN};
+Point(106) = {0.713778, 0.0000, -0.021289, msN};
+Point(107) = {0.735698, 0.0000, -0.018232, msN};
+Point(108) = {0.757051, 0.0000, -0.015357, msN};
+Point(109) = {0.777785, 0.0000, -0.012690, msN};
+Point(110) = {0.797850, 0.0000, -0.010244, msN};
+Point(111) = {0.817197, 0.0000, -0.008027, msN};
+Point(112) = {0.835779, 0.0000, -0.006048, msN};
+Point(113) = {0.853553, 0.0000, -0.004314, msN};
+Point(114) = {0.870476, 0.0000, -0.002829, msN};
+Point(115) = {0.886505, 0.0000, -0.001592, msN};
+Point(116) = {0.901604, 0.0000, -0.000600, msN};
+Point(117) = {0.915735, 0.0000, 0.000157, msN};
+Point(118) = {0.928864, 0.0000, 0.000694, msN};
+Point(119) = {0.940961, 0.0000, 0.001033, msN};
+Point(120) = {0.951995, 0.0000, 0.001197, msN};
+Point(121) = {0.961940, 0.0000, 0.001212, msN};
+Point(122) = {0.970772, 0.0000, 0.001112, msN};
+Point(123) = {0.978470, 0.0000, 0.000935, msN};
+Point(124) = {0.985016, 0.0000, 0.000719, msN};
+Point(125) = {0.990393, 0.0000, 0.000497, msN};
+Point(126) = {0.994588, 0.0000, 0.000296, msN};
+Point(127) = {0.997592, 0.0000, 0.000137, msN};
+Point(128) = {0.999398, 0.0000, 0.000035, msN};
+
+// Airfoil 2
+Point(301) = {1.000000, spn, 0.000000, msN};
+Point(302) = {0.999398, spn, 0.000128, msN};
+Point(303) = {0.997592, spn, 0.000510, msN};
+Point(304) = {0.994588, spn, 0.001137, msN};
+Point(305) = {0.990393, spn, 0.002001, msN};
+Point(306) = {0.985016, spn, 0.003092, msN};
+Point(307) = {0.978470, spn, 0.004401, msN};
+Point(308) = {0.970772, spn, 0.005915, msN};
+Point(309) = {0.961940, spn, 0.007622, msN};
+Point(310) = {0.951995, spn, 0.009508, msN};
+Point(311) = {0.940961, spn, 0.011562, msN};
+Point(312) = {0.928864, spn, 0.013769, msN};
+Point(313) = {0.915735, spn, 0.016113, msN};
+Point(314) = {0.901604, spn, 0.018580, msN};
+Point(315) = {0.886505, spn, 0.021153, msN};
+Point(316) = {0.870476, spn, 0.023817, msN};
+Point(317) = {0.853553, spn, 0.026554, msN};
+Point(318) = {0.835779, spn, 0.029347, msN};
+Point(319) = {0.817197, spn, 0.032176, msN};
+Point(320) = {0.797850, spn, 0.035017, msN};
+Point(321) = {0.777785, spn, 0.037847, msN};
+Point(322) = {0.757051, spn, 0.040641, msN};
+Point(323) = {0.735698, spn, 0.043377, msN};
+Point(324) = {0.713778, spn, 0.046029, msN};
+Point(325) = {0.691342, spn, 0.048575, msN};
+Point(326) = {0.668445, spn, 0.050993, msN};
+Point(327) = {0.645142, spn, 0.053258, msN};
+Point(328) = {0.621490, spn, 0.055344, msN};
+Point(329) = {0.597545, spn, 0.057218, msN};
+Point(330) = {0.573365, spn, 0.058845, msN};
+Point(331) = {0.549009, spn, 0.060194, msN};
+Point(332) = {0.524534, spn, 0.061254, msN};
+Point(333) = {0.500000, spn, 0.062029, msN};
+Point(334) = {0.475466, spn, 0.062530, msN};
+Point(335) = {0.450991, spn, 0.062774, msN};
+Point(336) = {0.426635, spn, 0.062779, msN};
+Point(337) = {0.402455, spn, 0.062562, msN};
+Point(338) = {0.378510, spn, 0.062133, msN};
+Point(339) = {0.354858, spn, 0.061497, msN};
+Point(340) = {0.331555, spn, 0.060660, msN};
+Point(341) = {0.308658, spn, 0.059629, msN};
+Point(342) = {0.286222, spn, 0.058414, msN};
+Point(343) = {0.264302, spn, 0.057026, msN};
+Point(344) = {0.242949, spn, 0.055470, msN};
+Point(345) = {0.222215, spn, 0.053753, msN};
+Point(346) = {0.202150, spn, 0.051885, msN};
+Point(347) = {0.182803, spn, 0.049874, msN};
+Point(348) = {0.164221, spn, 0.047729, msN};
+Point(349) = {0.146447, spn, 0.045457, msN};
+Point(350) = {0.129524, spn, 0.043071, msN};
+Point(351) = {0.113495, spn, 0.040585, msN};
+Point(352) = {0.098396, spn, 0.038011, msN};
+Point(353) = {0.084265, spn, 0.035360, msN};
+Point(354) = {0.071136, spn, 0.032644, msN};
+Point(355) = {0.059039, spn, 0.029874, msN};
+Point(356) = {0.048005, spn, 0.027062, msN};
+Point(357) = {0.038060, spn, 0.024219, msN};
+Point(358) = {0.029228, spn, 0.021348, msN};
+Point(359) = {0.021530, spn, 0.018441, msN};
+Point(360) = {0.014984, spn, 0.015489, msN};
+Point(361) = {0.009607, spn, 0.012480, msN};
+Point(362) = {0.005412, spn, 0.009416, msN};
+Point(363) = {0.002408, spn, 0.006306, msN};
+Point(364) = {0.000602, spn, 0.003165, msN};
+Point(365) = {0.000000, spn, 0.000000, msN};
+Point(366) = {0.000602, spn, -0.003160, msN};
+Point(367) = {0.002408, spn, -0.006308, msN};
+Point(368) = {0.005412, spn, -0.009443, msN};
+Point(369) = {0.009607, spn, -0.012559, msN};
+Point(370) = {0.014984, spn, -0.015649, msN};
+Point(371) = {0.021530, spn, -0.018707, msN};
+Point(372) = {0.029228, spn, -0.021722, msN};
+Point(373) = {0.038060, spn, -0.024685, msN};
+Point(374) = {0.048005, spn, -0.027586, msN};
+Point(375) = {0.059039, spn, -0.030416, msN};
+Point(376) = {0.071136, spn, -0.033170, msN};
+Point(377) = {0.084265, spn, -0.035843, msN};
+Point(378) = {0.098396, spn, -0.038431, msN};
+Point(379) = {0.113495, spn, -0.040929, msN};
+Point(380) = {0.129524, spn, -0.043326, msN};
+Point(381) = {0.146447, spn, -0.045610, msN};
+Point(382) = {0.164221, spn, -0.047773, msN};
+Point(383) = {0.182803, spn, -0.049805, msN};
+Point(384) = {0.202150, spn, -0.051694, msN};
+Point(385) = {0.222215, spn, -0.053427, msN};
+Point(386) = {0.242949, spn, -0.054994, msN};
+Point(387) = {0.264302, spn, -0.056376, msN};
+Point(388) = {0.286222, spn, -0.057547, msN};
+Point(389) = {0.308658, spn, -0.058459, msN};
+Point(390) = {0.331555, spn, -0.059046, msN};
+Point(391) = {0.354858, spn, -0.059236, msN};
+Point(392) = {0.378510, spn, -0.058974, msN};
+Point(393) = {0.402455, spn, -0.058224, msN};
+Point(394) = {0.426635, spn, -0.056979, msN};
+Point(395) = {0.450991, spn, -0.055257, msN};
+Point(396) = {0.475466, spn, -0.053099, msN};
+Point(397) = {0.500000, spn, -0.050563, msN};
+Point(398) = {0.524534, spn, -0.047719, msN};
+Point(399) = {0.549009, spn, -0.044642, msN};
+Point(400) = {0.573365, spn, -0.041397, msN};
+Point(401) = {0.597545, spn, -0.038043, msN};
+Point(402) = {0.621490, spn, -0.034631, msN};
+Point(403) = {0.645142, spn, -0.031207, msN};
+Point(404) = {0.668445, spn, -0.027814, msN};
+Point(405) = {0.691342, spn, -0.024495, msN};
+Point(406) = {0.713778, spn, -0.021289, msN};
+Point(407) = {0.735698, spn, -0.018232, msN};
+Point(408) = {0.757051, spn, -0.015357, msN};
+Point(409) = {0.777785, spn, -0.012690, msN};
+Point(410) = {0.797850, spn, -0.010244, msN};
+Point(411) = {0.817197, spn, -0.008027, msN};
+Point(412) = {0.835779, spn, -0.006048, msN};
+Point(413) = {0.853553, spn, -0.004314, msN};
+Point(414) = {0.870476, spn, -0.002829, msN};
+Point(415) = {0.886505, spn, -0.001592, msN};
+Point(416) = {0.901604, spn, -0.000600, msN};
+Point(417) = {0.915735, spn, 0.000157, msN};
+Point(418) = {0.928864, spn, 0.000694, msN};
+Point(419) = {0.940961, spn, 0.001033, msN};
+Point(420) = {0.951995, spn, 0.001197, msN};
+Point(421) = {0.961940, spn, 0.001212, msN};
+Point(422) = {0.970772, spn, 0.001112, msN};
+Point(423) = {0.978470, spn, 0.000935, msN};
+Point(424) = {0.985016, spn, 0.000719, msN};
+Point(425) = {0.990393, spn, 0.000497, msN};
+Point(426) = {0.994588, spn, 0.000296, msN};
+Point(427) = {0.997592, spn, 0.000137, msN};
+Point(428) = {0.999398, spn, 0.000035, msN};
+
+// Box
+Point(10001) = {1+lgt/2, 0, 0, msF};
+Point(10002) = {1+lgt/2, 0, hgt/2, msF};
+Point(10003) = {-lgt/2, 0, hgt/2, msF};
+Point(10004) = {-lgt/2, 0, -hgt/2, msF};
+Point(10005) = {1+lgt/2, 0, -hgt/2, msF};
+Point(10011) = {1+lgt/2, wdt, hgt/2, msF};
+Point(10012) = {-lgt/2, wdt, hgt/2, msF};
+Point(10013) = {-lgt/2, wdt, -hgt/2, msF};
+Point(10014) = {1+lgt/2, wdt, -hgt/2, msF};
+
+// Wake
+Point(10021) = {1+lgt/2, spn, 0, msF};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1:14};
+Spline(2) = {14:46};
+Spline(3) = {46:65};
+Spline(4) = {65:84};
+Spline(5) = {84:116};
+Spline(6) = {116:128,1};
+
+// Airfoil 2:
+Spline(7) = {301:314};
+Spline(8) = {314:346};
+Spline(9) = {346:365};
+Spline(10) = {365:384};
+Spline(11) = {384:416};
+Spline(12) = {416:428,301};
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,301};
+Line(62) = {14,314};
+Line(63) = {46,346};
+Line(64) = {65,365};
+Line(65) = {84,384};
+Line(66) = {116,416};
+
+// Box
+Line(201) = {10001,10002};
+Line(202) = {10002,10003};
+Line(203) = {10003,10004};
+Line(204) = {10004,10005};
+Line(205) = {10005,10001};
+Line(211) = {10011,10012};
+Line(212) = {10012,10013};
+Line(213) = {10013,10014};
+Line(214) = {10014,10011};
+Line(221) = {1, 10001};
+Line(222) = {301, 10021};
+Line(231) = {10002,10011};
+Line(232) = {10003,10012};
+Line(233) = {10004,10013};
+Line(234) = {10005,10014};
+
+// Wake
+Line(241) = {10001, 10021};
+
+/// Line loops & Surfaces
+
+// Wing 1:
+Line Loop(1) = {1,62,-7,-61};
+Line Loop(2) = {2,63,-8,-62};
+Line Loop(3) = {3,64,-9,-63};
+Line Loop(4) = {4,65,-10,-64};
+Line Loop(5) = {5,66,-11,-65};
+Line Loop(6) = {6,61,-12,-66};
+Surface(1) = {-1};
+Surface(2) = {-2};
+Surface(3) = {-3};
+Surface(4) = {-4};
+Surface(5) = {-5};
+Surface(6) = {-6};
+
+// Wingtip:
+Line Loop(11) = {7:12};
+Plane Surface(11) = {-11};
+
+// Symmetry
+Line Loop(21) = {201, 202, 203, 204, 205};
+Line Loop(22) = {1, 2, 3, 4, 5, 6};
+Plane Surface(21) = {-21, 22};
+
+// Downsteam
+Line Loop(31) = {201, 231, -214, -234, 205};
+Plane Surface(31) = {31};
+
+// Farfield
+Line Loop(41) = {233, -212, -232, 203};
+Plane Surface(41) = {41};
+Line Loop(42) = {204, 234, -213, -233};
+Plane Surface(42) = {42};
+Line Loop(43) = {202, 232, -211, -231};
+Plane Surface(43) = {43};
+Line Loop(44) = {213, 214, 211, 212};
+Plane Surface(44) = {44};
+
+// Wake
+Line Loop(51) = {221, 241, -222, -61};
+Surface(51) = {51};
+
+/// Volume
+Surface Loop(1) = {1,2,3,4,5,6,11,21,31,41,42,43,44};
+Volume(1) = {1};
+
+/// Embbeded
+Line{221} In Surface{21};
+Line{241} In Surface{31};
+Surface{51} In Volume{1};
+
+//// MESHING ALGORITHM
+
+/// 2D:
+
+///Wing, farfield and symmetry plane:
+MeshAlgorithm Surface{1,2,3,4,5,6,11,21,31,41,42,43,44,51} = 5;
+
+/// 3D:
+
+Mesh.Algorithm3D = 2;
+//Mesh.OptimizeNetgen = 1;
+Mesh.Optimize = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {222};
+Physical Line("te") = {61};
+
+// Internal Field:
+Physical Volume("field") = {1};
+
+// Wing:
+Physical Surface("wing") = {1,2,3,11};
+Physical Surface("wing_") = {4,5,6};
+
+// Symmetry:
+Physical Surface("symmetry") = {21};
+
+// Downstream:
+Physical Surface("downstream") = {31};
+
+// Farfield:
+Physical Surface("upstream") = {41};
+Physical Surface("farfield") = {42,43,44};
+
+// Wake:
+Physical Surface("wake") = {51};
diff --git a/dart/models/wbht.geo b/dart/models/wbht.geo
index 0e79a8e..9f252d8 100644
--- a/dart/models/wbht.geo
+++ b/dart/models/wbht.geo
@@ -950,9 +950,9 @@ Line(3026) = {3049, 3149};
/// Wake
// Wing
-Point(4001) = {25, 0, -2.26};
-Point(4002) = {25, 6.1, -0.66};
-Point(4003) = {25, 7, -0.41};
+Point(4001) = {25, 0, -0.534};
+Point(4002) = {25, 6.1, 0.01642};
+Point(4003) = {25, 7, 0.2577};
Line(4001) = {2001, 4001};
Line(4002) = {2301, 4002};
@@ -1213,10 +1213,12 @@ Mesh.SmoothNormals = 1;
/* Physical */
// Tips:
-Physical Line("teTipW") = {4092,4091,4094,2041,4003};
+Physical Line("teW") = {4094,2041};
Physical Line("wakeTipW") = {4003};
-Physical Line("teTipT") = {4103,4102,4101,4100,4105,4106};
+Physical Line("wakeExW") = {4092, 4091, 4003};
+Physical Line("teT") = {4105};
Physical Line("wakeTipT") = {4106};
+Physical Line("wakeExT") = {4102,4101,4100, 4103, 4106};
// Symmetry:
Physical Surface("symmetry") = {1};
diff --git a/dart/models/wht.geo b/dart/models/wht.geo
index 53a64ea..4fb94de 100644
--- a/dart/models/wht.geo
+++ b/dart/models/wht.geo
@@ -30,7 +30,7 @@ zOfst = 0.56;
// Domain
Point(1) = {-10, 0, -10};
Point(2) = {17, 0, -10};
-Point(4001) = {17, 0, -1.39};
+Point(4001) = {17, 0, -0.03763};
Point(4011) = {17, 0, 1.06};
Point(3) = {17, 0, 10};
Point(4) = {-10, 0, 10};
@@ -833,8 +833,8 @@ Line(3026) = {3049, 3149};
/// Wake
// Wing
-Point(4002) = {17, 6.1, 0.09};
-Point(4003) = {17, 7, 0.29};
+Point(4002) = {17, 6.1, 0.5764};
+Point(4003) = {17, 7, 0.8177};
Line(4001) = {2001, 4001};
Line(4002) = {2301, 4002};
@@ -988,9 +988,9 @@ Mesh.SmoothNormals = 1;
/* Physical */
// Tips:
-Physical Line("teTipW") = {2031, 2041, 4003};
+Physical Line("teW") = {2031, 2041};
Physical Line("wakeTipW") = {4003};
-Physical Line("teTipT") = {3021, 4012};
+Physical Line("teT") = {3021};
Physical Line("wakeTipT") = {4012};
// Symmetry:
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 00e62aa..23fe93f 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -117,8 +117,8 @@ void Adjoint::run()
dCf_X[1][i] -= dCd_XU[i]; // dCd/dX - dRu/dX^T * lambdaCdU
}
// compute mesh adjoint (totals)
- std::vector<double> lamClMsh = this->solveMesh(dCf_X[0]); // lambdaClX = dRx/dX^T * (dCl/dX - dRu/dX^T * lambdaClU)
- std::vector<double> lamCdMsh = this->solveMesh(dCf_X[1]); // lambdaCdX = dRx/dX^T * (dCd/dX - dRu/dX^T * lambdaCdU)
+ std::vector<double> lamClMsh = this->solveMesh(dCf_X[0]); // lambdaClX = dRx/dX^-T * (dCl/dX - dRu/dX^T * lambdaClU)
+ std::vector<double> lamCdMsh = this->solveMesh(dCf_X[1]); // lambdaCdX = dRx/dX^-T * (dCd/dX - dRu/dX^T * lambdaCdU)
// store solution
for (auto n : sol->pbl->msh->nodes)
{
diff --git a/dart/src/wKuttaElement.cpp b/dart/src/wKuttaElement.cpp
index c66a8de..360efa5 100644
--- a/dart/src/wKuttaElement.cpp
+++ b/dart/src/wKuttaElement.cpp
@@ -28,26 +28,57 @@ KuttaElement::KuttaElement(size_t _no, Element *&_surE, Element *&_volE,
// Sanity checks
if (surE->type() == ELTYPE::LINE2)
nDim = 2;
+ else if (surE->type() == ELTYPE::TRI3)
+ nDim = 3;
else
{
std::stringstream err;
err << "KuttaElement only implemented for surface elements of type "
- << ELTYPE::LINE2 << " (" << surE->type() << " was given)!\n";
+ << ELTYPE::LINE2 << " or " << ELTYPE::TRI3
+ << " (" << surE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
- if (volE->type() == ELTYPE::TRI3)
+ if (volE->type() == ELTYPE::TRI3 || volE->type() == ELTYPE::TETRA4)
nCol = volE->nodes.size();
else
{
std::stringstream err;
err << "KuttaElement only implemented for volume elements of type "
- << ELTYPE::TRI3 << " (" << volE->type() << " was given)!\n";
+ << ELTYPE::TRI3 << " or " << ELTYPE::TETRA4
+ << " (" << volE->type() << " was given)!\n";
throw std::runtime_error(err.str());
}
// Count TE nodes
+ mapNodes();
nRow = teN.size();
}
+/**
+ * @brief Map the wake nodes to their local volume node indices
+ */
+void KuttaElement::mapNodes()
+{
+ // Find local volume index of contributing surface nodes
+ for (auto p : teN)
+ {
+ for (size_t i = 0; i < volE->nodes.size(); ++i)
+ {
+ if (surE->nodes[p.first] == volE->nodes[i])
+ {
+ vsMap[p.first] = i;
+ break;
+ }
+ }
+ }
+ // Check that all contributing surface nodes have been mapped
+ if (vsMap.size() != teN.size())
+ {
+ std::stringstream err;
+ err << "dart::KuttaElement: could not map some contributing nodes from surface element " << *surE << " to " << *volE << " (number of mapped nodes: " << vsMap.size() << ")!\n";
+ throw std::runtime_error(err.str());
+ }
+}
+
void KuttaElement::write(std::ostream &out) const
{
out << "KuttaElement #" << no
diff --git a/dart/src/wKuttaElement.h b/dart/src/wKuttaElement.h
index b48c48c..f3fe8d8 100644
--- a/dart/src/wKuttaElement.h
+++ b/dart/src/wKuttaElement.h
@@ -21,6 +21,7 @@
#include "wObject.h"
#include "wElement.h"
#include <vector>
+#include <map>
namespace dart
{
@@ -29,16 +30,12 @@ namespace dart
* @brief Super-element to apply Kutta condition
*
* Made up of one surface element at the trailing edge and the connected volume
- * element. Only Line2 (surface) and Tri3 (volume) elements are curently
- * implemented.
- * Since Tri3 are used, the evaluation of the volume integrals is performed at
- * dummy Gauss point #0. If volume elements with non-constant shape function
- * derivatives are used, those derivates should be recomputed at the true
- * surface element Gauss points.
- * Kutta works well in 2D, but fails in 3D unless the geometry follows strict
- * restrictions (thin and sharp TE, simple planform, same mesh elements on TE
- * suction and pressure side, ...). Kutta contribution is therefore disabled for
- * 3D cases.
+ * element. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are
+ * currently implemented.
+ * Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is
+ * performed at dummy Gauss point #0. If volume elements with non-constant shape
+ * function derivatives are used, those derivates should be recomputed at the
+ * true surface element Gauss points.
* @authors Adrien Crovato
*/
class DART_API KuttaElement : public fwk::wObject
@@ -54,12 +51,16 @@ public:
size_t nDim; ///< dimension of volume element
std::vector<std::pair<size_t, tbox::Node *>> teN; ///< Map of local node indices and trailing edge nodes
+ std::map<size_t, size_t> vsMap; ///< map between surface node and volume node indices
int sign; ///< +1 if upper element, -1 if lower element
KuttaElement(size_t _no, tbox::Element *&_surE, tbox::Element *&_volE, std::vector<std::pair<size_t, tbox::Node *>> &_teN, int _sign);
virtual void write(std::ostream &out) const override;
#endif
+
+private:
+ void mapNodes();
};
} // namespace dart
diff --git a/dart/src/wKuttaResidual.cpp b/dart/src/wKuttaResidual.cpp
index cbca6a6..6d0b9e3 100644
--- a/dart/src/wKuttaResidual.cpp
+++ b/dart/src/wKuttaResidual.cpp
@@ -27,18 +27,30 @@ using namespace dart;
/**
* @brief Build the residual matrix for a fixed-point iteration, on one Kutta element
*
- * A = \int psi * v_u*grad_phi dS
+ * A = 1/S * \int (psi + grad_psi*v_inf) * v*grad_phi dS
+ * NB: v_inf = [1, 0, 0]
*/
Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<double> const &phi)
{
- // Get shape functions, Jacobian and Gauss points
+ // Get shape functions and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
+ Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
+ // Get Jacobian
+ Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
// Get size
size_t nRow = ke.nRow;
+ size_t nDim = ke.nDim;
+ // Stabilization
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t j = 0; j < nDim; ++j)
+ dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first));
+ Eigen::VectorXd dSfp(nRow);
+ dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
// Velocity
- Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getJinv(0) * ke.volE->getVCache().getDsf(0);
+ Eigen::RowVectorXd Aup = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).transpose() * iJ * dSfV;
// Build
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
@@ -48,7 +60,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
- sfp(i) = sf(ke.teN[i].first);
+ sfp(i) = sf(ke.teN[i].first) + dSfp(i);
// wake contribution
A += ke.sign * sfp * Aup * sGauss.getW(k) * ke.surE->getDetJ(k);
}
@@ -58,19 +70,30 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
/**
* @brief Build the residual vector, on one Kutta element
*
- * b = \int psi * (grad_phi)^2 dS
+ * b = 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS
+ * NB: v_inf = [1, 0, 0]
*/
Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> const &phi)
{
- // Get shape functions, Jacobian and Gauss points
+ // Get shape functions and Gauss points
Cache &sCache = ke.surE->getVCache();
Gauss &sGauss = sCache.getVGauss();
+ Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
+ // Get Jacobian
+ Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
// Get size
size_t nRow = ke.nRow;
+ size_t nDim = ke.nDim;
+ // Stabilization
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t j = 0; j < nDim; ++j)
+ dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first));
+ Eigen::VectorXd dSfp(nRow);
+ dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
// Velocity
- Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
- double dPhi2 = dPhi.squaredNorm();
+ double dPhi2 = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).squaredNorm();
// Build
Eigen::VectorXd b = Eigen::VectorXd::Zero(nRow);
@@ -80,7 +103,7 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
- sfp(i) = sf(ke.teN[i].first);
+ sfp(i) = sf(ke.teN[i].first) + dSfp(i);
// wake contribution
b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * ke.surE->getDetJ(k);
}
@@ -90,7 +113,8 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
/**
* @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one Kutta element
*
- * A = \int psi * 2*(grad_phi_u)*dgrad_phi_u dS
+ * A = 1/S * \int (psi + grad_psi*v_inf) * 2*(grad_phi_u)*dgrad_phi_u dS
+ * NB: v_inf = [1, 0, 0]
*/
Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::vector<double> const &phi)
{
@@ -101,27 +125,62 @@ Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::ve
/**
* @brief Build the gradient of the residual vector with respect to the nodes, on one Kutta element
- *
- * A = d( \int psi * (grad_phi)^2 dS )
- * = \int psi * d(grad_phi)^2 dS
- * + \int psi * (grad_phi)^2 ddS
+ *
+ * A = d( 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS )
+ * = d1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS
+ * + 1/S * \int dgrad_psi*v_inf * (grad_phi)^2 dS
+ * + 1/S * \int (psi + grad_psi*v_inf) * d(grad_phi)^2 dS
+ * + 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 ddS
+ * NB: v_inf = [1, 0, 0]
*/
std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(KuttaElement const &ke, std::vector<double> const &phi)
{
// Get shape functions and Gauss points
Cache &surCache = ke.surE->getVCache();
Gauss &surGauss = surCache.getVGauss();
+ Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0);
// Get Jacobian
Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
std::vector<Eigen::MatrixXd> const &dJ = ke.volE->getGradJ(0);
+ // Get surface and gradient
+ double surf = ke.surE->getVol();
+ std::vector<double> const &dSurf = ke.surE->getGradVol();
// Get size
size_t nRow = ke.nRow;
size_t nDim = ke.nDim;
size_t nCol = ke.nCol;
size_t nSur = ke.surE->nodes.size();
- // velocity
+ // Stabilization term
+ Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
+ double sqSurf = sqrt(surf);
+ Eigen::MatrixXd dSf(nRow, nDim);
+ for (size_t i = 0; i < nRow; ++i)
+ for (size_t m = 0; m < nDim; ++m)
+ dSf(i, m) = dSfV(m, ke.vsMap.at(ke.teN[i].first));
+ Eigen::VectorXd gradSfp = 0.5 * sqSurf * dSf * iJ.transpose() * vi;
+ // Gradient of stabilization term (volume and surface)
+ std::vector<Eigen::VectorXd> dGradSfV(nDim * nCol);
+ for (size_t i = 0; i < nCol; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi
+ }
+ }
+ std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur);
+ for (size_t i = 0; i < nSur; ++i)
+ {
+ for (size_t m = 0; m < nDim; ++m)
+ {
+ size_t idx = i * nDim + m;
+ dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJ.transpose() * vi; // dS * grad_psi
+ }
+ }
+ // Velocity and normalization
Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0);
+ double nrm = 1 / surf;
// Build
Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(nRow, nDim * nCol);
@@ -136,7 +195,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
Eigen::VectorXd const &sf = surCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
- sfp(i) = sf(ke.teN[i].first);
+ sfp(i) = sf(ke.teN[i].first) + gradSfp(i);
// Volume
for (size_t i = 0; i < nCol; ++i)
@@ -145,7 +204,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
{
size_t idx = i * nDim + m;
// gradient of velocity jump
- Av.col(idx) += ke.sign * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // psi * dgrad_phi^2
+ Av.col(idx) += ke.sign * nrm * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // 1/S * (psi + grad_psi) * dgrad_phi^2 * detJ
+ // gradient of stabilization term (volume)
+ Av.col(idx) += ke.sign * nrm * w * dGradSfV[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ
}
}
// Surface
@@ -154,8 +215,12 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
for (size_t m = 0; m < nDim; ++m)
{
size_t idx = i * nDim + m;
+ // gradient of normalization
+ As.col(idx) -= ke.sign * nrm * nrm * dSurf[idx] * w * sfp * dPhi.dot(dPhi) * detJ; // d1/S * (psi + grad_psi) * grad_phi^2 * detJ
+ // gradient of stabilization term (surface)
+ As.col(idx) += ke.sign * nrm * w * dGradSfS[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ
// gradient of jacobian determinant
- As.col(idx) += ke.sign * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // psi * grad_phi^2 * ddetJ
+ As.col(idx) += ke.sign * nrm * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // 1/S * (psi + grad_psi) * grad_phi^2 * ddetJ
}
}
}
diff --git a/dart/src/wProblem.cpp b/dart/src/wProblem.cpp
index 73d6f91..06d219d 100644
--- a/dart/src/wProblem.cpp
+++ b/dart/src/wProblem.cpp
@@ -194,63 +194,42 @@ void Problem::check() const
// Sanity checks
if (fluid == NULL)
throw std::runtime_error("No Fluid provided!\n");
+ if (bodies.empty())
+ throw std::runtime_error("No Body provided!\n");
if (iIC == NULL)
throw std::runtime_error("No Initial Condition provided!\n");
if (fBCs.empty())
throw std::runtime_error("No Freestream B.C. provided!\n");
if (wBCs.empty())
throw std::runtime_error("No Wake B.C. provided!\n");
+ if (kSCs.empty())
+ throw std::runtime_error("No Kutta condition provided!\n");
// Three-dimension problem
if (nDim == 3)
{
// Volume element type
for (auto e : fluid->tag->elems)
- {
if (e->type() != ELTYPE::TETRA4)
- {
- std::stringstream err;
- err << "3D solver is only implemented for volume elements of type "
- << ELTYPE::TETRA4 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "3", "volume", ELTYPE::TETRA4);
+ // Surface element type (Body)
+ for (auto surf : bodies)
+ for (auto e : surf->groups[0]->tag->elems)
+ if (e->type() != ELTYPE::TRI3)
+ throwUnsupElType(e, "3", "surface", ELTYPE::TRI3);
// Surface element type (Freestream B.C.)
for (auto surf : fBCs)
- {
for (auto e : surf->tag->elems)
- {
if (e->type() != ELTYPE::TRI3)
- {
- std::stringstream err;
- err << "3D solver is only implemented for surface elements of type "
- << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
- }
+ throwUnsupElType(e, "3", "surface", ELTYPE::TRI3);
// Surface element type (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
- {
if (e->type() != ELTYPE::TRI3)
- {
- std::stringstream err;
- err << "3D solver is only implemented for surface elements of type "
- << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "3", "surface", ELTYPE::TRI3);
for (auto e : surf->groups[1]->tag->elems)
- {
if (e->type() != ELTYPE::TRI3)
- {
- std::stringstream err;
- err << "3D solver is only implemented for surface elements of type "
- << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "3", "surface", ELTYPE::TRI3);
}
// Blowing B.C.
if (!bBCs.empty())
@@ -261,52 +240,27 @@ void Problem::check() const
{
// Volume element type
for (auto e : fluid->tag->elems)
- {
if (e->type() != ELTYPE::TRI3)
- {
- std::stringstream err;
- err << "2D solver is only implemented for volume elements of type "
- << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "2", "volume", ELTYPE::TRI3);
+ // Surface element type (Body)
+ for (auto surf : bodies)
+ for (auto e : surf->groups[0]->tag->elems)
+ if (e->type() != ELTYPE::LINE2)
+ throwUnsupElType(e, "2", "surface", ELTYPE::LINE2);
// Surface element type (Freestream B.C.)
for (auto surf : fBCs)
- {
for (auto e : surf->tag->elems)
- {
if (e->type() != ELTYPE::LINE2)
- {
- std::stringstream err;
- err << "2D solver is only implemented for surface elements of type "
- << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
- }
+ throwUnsupElType(e, "2", "surface", ELTYPE::LINE2);
// Surface element type (Wake B.C.)
for (auto surf : wBCs)
{
for (auto e : surf->groups[0]->tag->elems)
- {
if (e->type() != ELTYPE::LINE2)
- {
- std::stringstream err;
- err << "2D solver is only implemented for surface elements of type "
- << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "2", "surface", ELTYPE::LINE2);
for (auto e : surf->groups[1]->tag->elems)
- {
if (e->type() != ELTYPE::LINE2)
- {
- std::stringstream err;
- err << "2D solver is only implemented for surface elements of type "
- << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n";
- throw std::runtime_error(err.str());
- }
- }
+ throwUnsupElType(e, "2", "surface", ELTYPE::LINE2);
}
}
else
@@ -318,6 +272,15 @@ void Problem::check() const
}
}
+void Problem::throwUnsupElType(Element const *e, std::string const &pdim, std::string const &edim, ELTYPE type) const
+{
+ std::stringstream err;
+ err << pdim << "D solver is only implemented for "
+ << edim << " elements of type " << type
+ << " (" << e->type() << " was given)!\n";
+ throw std::runtime_error(err.str());
+}
+
void Problem::write(std::ostream &out) const
{
out << "dart::Problem parameters"
diff --git a/dart/src/wProblem.h b/dart/src/wProblem.h
index 3ba0200..d59911e 100644
--- a/dart/src/wProblem.h
+++ b/dart/src/wProblem.h
@@ -19,6 +19,7 @@
#include "dart.h"
#include "wObject.h"
+#include "wElement.h"
#include "wF1Ct.h"
#include <memory>
#include <iostream>
@@ -87,6 +88,9 @@ public:
void initGradElems();
virtual void write(std::ostream &out) const override;
#endif
+
+private:
+ void throwUnsupElType(tbox::Element const *e, std::string const &pdim, std::string const &edim, tbox::ELTYPE type) const;
};
} // namespace dart
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 7f69f84..3b11893 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -21,6 +21,7 @@
#include "wAssign.h"
#include "wBody.h"
#include "wWake.h"
+#include "wKutta.h"
#include "wLoadFunctional.h"
#include "wMshData.h"
@@ -110,6 +111,8 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
std::cout << "Dirichlet on " << *bnd->tag << "\n";
for (auto bnd : pbl->wBCs)
std::cout << "Wake on " << *bnd->groups[0]->tag << "\n";
+ for (auto bnd : pbl->kSCs)
+ std::cout << "Kutta on " << *bnd->groups[0]->tag << "\n";
for (auto bnd : pbl->bodies)
std::cout << "Body on " << *bnd->groups[0]->tag << "\n";
std::cout << "--- Freestream conditions ---\n"
diff --git a/dart/src/wWakeElement.h b/dart/src/wWakeElement.h
index 0cf0a30..43ead69 100644
--- a/dart/src/wWakeElement.h
+++ b/dart/src/wWakeElement.h
@@ -30,8 +30,8 @@ namespace dart
* @brief Super-element to apply wake boundary conditions
*
* Made up of two surface elements on the wake and the connected volume
- * elements. Only Line2 (surface) and Tri3 (volume) elements are curently
- * implemented.
+ * elements. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are
+ * currently implemented.
* Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is
* performed at dummy Gauss point #0. If volume elements with non-constant shape
* function derivatives are used, those derivates should be recomputed at the
diff --git a/dart/src/wWakeResidual.cpp b/dart/src/wWakeResidual.cpp
index 671eb02..4af3176 100644
--- a/dart/src/wWakeResidual.cpp
+++ b/dart/src/wWakeResidual.cpp
@@ -66,7 +66,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
- sfp(i) = dSfp(i) + sf(we.wkN[i].first);
+ sfp(i) = sf(we.wkN[i].first) + dSfp(i);
// wake contribution
Au += sfp * Aup * sGauss.getW(k) * we.surUpE->getDetJ(k);
Al += sfp * Alw * sGauss.getW(k) * we.surUpE->getDetJ(k);
@@ -99,10 +99,8 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Eigen::VectorXd dSfp(nRow);
dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim);
// Velocity squared
- Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
- Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
- double dPhi2u = dPhiU.squaredNorm();
- double dPhi2l = dPhiL.squaredNorm();
+ double dPhi2u = we.volUpE->computeGradient(phi, 0).squaredNorm();
+ double dPhi2l = we.volLwE->computeGradient(phi, 0).squaredNorm();
// Build
Eigen::VectorXd bu = Eigen::VectorXd::Zero(nRow);
@@ -113,7 +111,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
Eigen::VectorXd const &sf = sCache.getSf(k);
Eigen::VectorXd sfp(nRow);
for (size_t i = 0; i < nRow; ++i)
- sfp(i) = dSfp(i) + sf(we.wkN[i].first);
+ sfp(i) = sf(we.wkN[i].first) + dSfp(i);
// wake contribution
bu += sfp * dPhi2u * sGauss.getW(k) * we.surUpE->getDetJ(k);
bl += sfp * dPhi2l * sGauss.getW(k) * we.surUpE->getDetJ(k);
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 3d3932a..6d35b84 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -55,7 +55,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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')
tms['pre'].stop()
# solve problem
@@ -114,17 +114,17 @@ def main():
if M_inf == 0. and alpha == 0*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3))
- tests.add(CTest('dCl_dMsh', dClX, 55.439, 1e-3))
- tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-3))
- tests.add(CTest('dCl_dAoA (msh)', dClAoA, 6.9, 1e-2))
+ tests.add(CTest('dCl_dMsh', dClX, 55.1, 1e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-2))
+ tests.add(CTest('dCl_dAoA (msh)', dClAoA, 7.1, 1e-2)) # FD 7.1364 (step = 1e-6)
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3))
elif M_inf == 0.7 and alpha == 2*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5)
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5)
- tests.add(CTest('dCl_dMsh', dClX, 82.67, 1e-3))
- tests.add(CTest('dCd_dMsh', dCdX, 0.908, 1e-3))
- tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2))
- tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2))
+ tests.add(CTest('dCl_dMsh', dClX, 83.6, 1e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.904, 1e-2))
+ tests.add(CTest('dCl_dAoA (msh)', dClAoA, 12.1, 1e-2)) # FD 12.1443 (step = 1e-6)
+ tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.129, 1e-2)) # FD 0.1289 (step = 1e-6)
else:
raise Exception('Test not defined for this flow')
tests.run()
diff --git a/dart/tests/bli.py b/dart/tests/bli.py
index fd2f186..9014bcc 100644
--- a/dart/tests/bli.py
+++ b/dart/tests/bli.py
@@ -57,7 +57,7 @@ def main():
M_inf = 0.5
c_ref = 1
dim = 2
- tol = 1e-4 # tolerance for the VII (usually one drag count)
+ tol = 1e-3 # tolerance for the VII (usually one drag count) #TODO tolerance has been decreased with new Kutta, hopefully will increase with Paul's dev
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
@@ -68,7 +68,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- 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)
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc = True, dbc=True)
tms['pre'].stop()
# solve viscous problem
diff --git a/dart/tests/cylinder.py b/dart/tests/cylinder.py
index 88dc36a..2d5f078 100644
--- a/dart/tests/cylinder.py
+++ b/dart/tests/cylinder.py
@@ -59,7 +59,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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', dbc=True)
tms['pre'].stop()
# solve problem
diff --git a/dart/tests/cylinder2D5.py b/dart/tests/cylinder2D5.py
deleted file mode 100644
index 107e6a4..0000000
--- a/dart/tests/cylinder2D5.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Compute the (non)linear flow around extruded 2D cylinder, with various B.C.
-# Adrien Crovato
-#
-# Test solver convergence and boundary conditions
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import math
-import dart.default as floD
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- alpha = 3*math.pi/180
- U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
- M_inf = 0.1
- dim = 3
-
- # define dimension and mesh size
- dmt = 1.0 # cylinder diameter
- lgt = 8.0 # channel length
- hgt = 6.0 # channel height
- wdt = 3.0 # channel width
- l_ref = dmt # reference length
- S_ref = dmt*wdt # reference area
- fms = 0.5 # farfield mesh size
- nms = 0.1 # nearfield mesh size
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms["msh"].start()
- pars = {'dmt' : dmt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/cylinder_2D5.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'])
- tms["msh"].stop()
-
- # set the problem and add fluid, initial/boundary conditions
- tms['pre'].start()
- 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
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
- solver = floD.newton(pbl)
- cnvrgd = solver.run()
- solver.save(gmshWriter)
- tms['solver'].stop()
-
- # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
- gradn = []
- for e in bnd.groups[0].tag.elems:
- gradx = 0
- grady = 0
- gradz = 0
- ux = e.nodes[1].pos[0] - e.nodes[0].pos[0]
- uy = e.nodes[1].pos[1] - e.nodes[0].pos[1]
- uz = e.nodes[1].pos[2] - e.nodes[0].pos[2]
- vx = e.nodes[2].pos[0] - e.nodes[0].pos[0]
- vy = e.nodes[2].pos[1] - e.nodes[0].pos[1]
- vz = e.nodes[2].pos[2] - e.nodes[0].pos[2]
- nx = uy*vz - uz*vy
- ny = uz*vx - ux*vz
- nz = ux*vy - uy*vx
- for n in e.nodes:
- gradx += solver.U[n.row][0] / e.nodes.size()
- grady += solver.U[n.row][1] / e.nodes.size()
- gradz += solver.U[n.row][2] / e.nodes.size()
- gradn_ = gradx * nx + grady * ny + gradz * nz
- gradn.append(gradn_)
- errNeu = sum(gradn) / len(gradn)
- # compute largest error in Dirichlet boundary condition
- errDir = 0.
- for n in dirichlet.nodes:
- errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2])
- if (abs(errCur) > errDir):
- errDir = abs(errCur)
- # compute largest error in Wake boundary condition
- errWak = 0.
- for n in wake.nodMap:
- errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row]
- if (abs(errCur) > errWak):
- errWak = abs(errCur)
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # visualize solution
- floD.initViewer(pbl)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- if (not cnvrgd):
- raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
- tests = CTests()
- tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
- tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2))
- tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12))
- tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1))
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/dart/tests/cylinder3.py b/dart/tests/cylinder3.py
deleted file mode 100644
index 869eca3..0000000
--- a/dart/tests/cylinder3.py
+++ /dev/null
@@ -1,137 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Compute the (non)linear flow around 3D cylinder, with various B.C.
-# Adrien Crovato
-#
-# Test solver convergence and boundary conditions
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import math
-import dart.utils as floU
-import dart.default as floD
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- alpha = 3*math.pi/180
- U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
- M_inf = 0.1
- dim = 3
-
- # define dimension and mesh size
- dmt = 1.0 # cylinder diameter
- spn = 2.0 # cylinder span
- tpt = 5e-2 # cylinder tip thickness (ratio wrt span)
- lgt = 8.0 # channel length
- hgt = 6.0 # channel height
- wdt = 3.0 # channel width
- l_ref = dmt # reference length
- S_ref = dmt*spn # reference area
- fms = 0.5 # farfield mesh size
- nms = 0.1 # nearfield mesh size
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'dmt' : dmt, 'spn' : spn, 'tpt' : tpt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/cylinder_3.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'], wktp = 'wakeTip')
- tms['msh'].stop()
-
- # set the problem and add fluid, initial/boundary conditions
- tms['pre'].start()
- 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
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
- solver = floD.newton(pbl)
- cnvrgd = solver.run()
- solver.save(gmshWriter)
- tms['solver'].stop()
-
- # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
- gradn = []
- for e in bnd.groups[0].tag.elems:
- gradx = 0
- grady = 0
- gradz = 0
- ux = e.nodes[1].pos[0] - e.nodes[0].pos[0]
- uy = e.nodes[1].pos[1] - e.nodes[0].pos[1]
- uz = e.nodes[1].pos[2] - e.nodes[0].pos[2]
- vx = e.nodes[2].pos[0] - e.nodes[0].pos[0]
- vy = e.nodes[2].pos[1] - e.nodes[0].pos[1]
- vz = e.nodes[2].pos[2] - e.nodes[0].pos[2]
- nx = uy*vz - uz*vy
- ny = uz*vx - ux*vz
- nz = ux*vy - uy*vx
- for n in e.nodes:
- gradx += solver.U[n.row][0] / e.nodes.size()
- grady += solver.U[n.row][1] / e.nodes.size()
- gradz += solver.U[n.row][2] / e.nodes.size()
- gradn_ = gradx * nx + grady * ny + gradz * nz
- gradn.append(gradn_)
- errNeu = sum(gradn) / len(gradn)
- # compute largest error in Dirichlet boundary condition
- errDir = 0.
- for n in dirichlet.nodes:
- errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2])
- if (abs(errCur) > errDir):
- errDir = abs(errCur)
- # compute largest error in Wake boundary condition
- errWak = 0.
- for n in wake.nodMap:
- errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row]
- if (abs(errCur) > errWak):
- errWak = abs(errCur)
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # visualize solution
- floD.initViewer(pbl)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- if (not cnvrgd):
- raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
- tests = CTests()
- tests.add(CTest('iteration count', solver.nIt, 4, 1, forceabs=True))
- tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2))
- tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12))
- tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1))
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/dart/tests/lift.py b/dart/tests/lift.py
index 54f58b3..1c8ee82 100644
--- a/dart/tests/lift.py
+++ b/dart/tests/lift.py
@@ -53,7 +53,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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')
tms['pre'].stop()
# solve problem
@@ -97,6 +97,7 @@ def main():
tests.add(CTest('iteration count', solver.nIt, 12, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.28, 5e-2))
tests.add(CTest('Cl', solver.Cl, 0.388, 5e-2))
+ tests.add(CTest('Cd', solver.Cd, 0.0013, 5e-4, forceabs=True))
tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
else:
raise Exception('Test not defined for this flow')
diff --git a/dart/tests/lift3.py b/dart/tests/lift3.py
deleted file mode 100644
index 67c2ad3..0000000
--- a/dart/tests/lift3.py
+++ /dev/null
@@ -1,101 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
-# Adrien Crovato
-#
-# Test the nonlinear shock-capturing capability and the Kutta condition
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import math
-import dart.utils as floU
-import dart.default as floD
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- alpha = 3*math.pi/180
- M_inf = 0.3
- dim = 3
-
- # define dimension and mesh size
- spn = 1.0 # wing span
- lgt = 3.0 # channel half length
- hgt = 3.0 # channel half height
- wdt = 3.0 # channel width
- S_ref = 1.*spn # reference area
- c_ref = 1 # reference chord
- fms = 1.0 # farfield mesh size
- nms = 0.02 # nearfield mesh size
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
- tms['msh'].stop()
-
- # set the problem and add fluid, initial/boundary conditions
- tms['pre'].start()
- 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
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
- solver = floD.newton(pbl)
- solver.run()
- solver.save(gmshWriter)
- tms['solver'].stop()
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # visualize solution
- floD.initViewer(pbl)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
- tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
- tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
- tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
- tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
- tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
- else:
- raise Exception('Test not defined for this flow')
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/dart/tests/meshDef3.py b/dart/tests/meshDef3.py
deleted file mode 100644
index d571ee6..0000000
--- a/dart/tests/meshDef3.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Compute the flow around a NACA 0012 wing on a deformed mesh
-# Adrien Crovato
-#
-# Test the mesh deformation process
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import numpy as np
-import dart.utils as floU
-import dart.default as floD
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- alpha0 = 0 # must be zero for this testfile
- M_inf = 0.3
- dim = 3
-
- # define dimension and mesh size
- spn = 1.0 # wing span
- lgt = 3.0 # channel half length
- hgt = 3.0 # channel half height
- wdt = 3.0 # channel width
- S_ref = 1.*spn # reference area
- c_ref = 1 # reference chord
- fms = 1.0 # farfield mesh size
- nms = 0.05 # nearfield mesh size
-
- # parameters for mesh deformation
- alfa = 3*np.pi/180
- xc = 0.25
- dz_max = 0.2*spn
-
- # mesh an airfoil
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
- tms['msh'].stop()
- # create mesh deformation handler
- mshDef = floD.morpher(msh, dim, ['wing'])
-
- # set the problem and add fluid, initial/boundary conditions
- tms['pre'].start()
- 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
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver0'].start()
- solver = floD.newton(pbl)
- solver.run()
- tms['solver0'].stop()
-
- # deform the mesh (dummy rotation + bending)
- print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET)
- tms['deform'].start()
- mshDef.savePos()
- rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]])
- vNods = {}
- for e in bnd.groups[0].tag.elems:
- for n in e.nodes:
- if n.no not in vNods:
- vNods[n.no] = n
- for no in vNods:
- nod = np.matrix([[vNods[no].pos[0]-xc], [vNods[no].pos[2]]])
- new = rMat * nod
- vNods[no].pos[0] = new[0,0]
- vNods[no].pos[2] = new[1,0] + dz_max/(spn*spn)*(vNods[no].pos[1]*vNods[no].pos[1])
- # deform the mesh
- mshDef.deform()
- gmshWriter.save(msh.name+"_def")
- tms['deform'].stop()
-
- # solve problem for alfa AOA
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver1'].start()
- solver.run()
- solver.save(gmshWriter)
- tms['solver1'].stop()
-
- # display results
- print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' case M alpha Cl Cd Cm')
- print(' true {0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, solver.Cm))
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- if alfa == 3*np.pi/180 and M_inf == 0.3 and spn == 1.0:
- tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
- tests.add(CTest('CL', solver.Cl, 0.135, 5e-1))
- tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
- tests.add(CTest('CS', solver.Cs, -0.011, 5e-2))
- tests.add(CTest('CM', solver.Cm, 0.0061, 1e-2))
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/dart/tests/meshDef.py b/dart/tests/morpher.py
similarity index 96%
rename from dart/tests/meshDef.py
rename to dart/tests/morpher.py
index c50e5af..c82f277 100644
--- a/dart/tests/meshDef.py
+++ b/dart/tests/morpher.py
@@ -55,7 +55,7 @@ def main():
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
- mshDef = floD.morpher(msh, dim, ['airfoil'])
+ morpher = floD.morpher(msh, dim, ['airfoil'])
# mesh the reference geometry
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01, 'angle' : alfa, 'xRot' : xc}
msh_ref, gmshWriter_ref = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
@@ -63,9 +63,9 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, 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')
# set the refrence problem and add fluid, initial/boundary conditions
- pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, 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')
tms['pre'].stop()
# solver problem for 0°AOA
@@ -78,7 +78,7 @@ def main():
# deform the mesh (dummy rotation + translation)
print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET)
tms['deform'].start()
- mshDef.savePos()
+ morpher.savePos()
rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]])
vNods = {}
for e in bnd.groups[0].tag.elems:
@@ -91,7 +91,7 @@ def main():
vNods[no].pos[0] = new[0,0] + dx + xc
vNods[no].pos[1] = new[1,0] + dy
# deform the mesh
- mshDef.deform()
+ morpher.deform()
gmshWriter.save("n0012_def")
tms['deform'].stop()
diff --git a/dart/tests/nonlift.py b/dart/tests/nonlift.py
index 0613660..2710fcb 100644
--- a/dart/tests/nonlift.py
+++ b/dart/tests/nonlift.py
@@ -54,7 +54,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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')
tms['pre'].stop()
# solve problem
@@ -93,6 +93,7 @@ def main():
elif M_inf == 0.8:
tests.add(CTest('iteration count', solver.nIt, 15, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
+ tests.add(CTest('Cd', solver.Cd, 0.0049, 5e-4, forceabs=True))
else:
raise Exception('Test not defined for this flow')
tests.add(CTest('Cl', solver.Cl, 0., 1e-2))
diff --git a/dart/tests/adjoint3.py b/dart/tests/rae_25.py
similarity index 62%
rename from dart/tests/adjoint3.py
rename to dart/tests/rae_25.py
index a9bb5a7..1647323 100644
--- a/dart/tests/adjoint3.py
+++ b/dart/tests/rae_25.py
@@ -2,13 +2,13 @@
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
-#
+#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
-#
+#
# http://www.apache.org/licenses/LICENSE-2.0
-#
+#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -16,22 +16,20 @@
# limitations under the License.
-## Compute adjoint solution of lifting (linear or nonlinear) potential flow around a NACA 0012
+## Compute the nonlinear flow around extruded 2D rae2822
# Adrien Crovato
#
-# Test the adjoint solver
+# Test transonic capabilities and kutta condition
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import numpy as np
-import dart.utils as floU
import dart.default as floD
-import tbox.utils as tboxU
import fwk
from fwk.testing import *
-from fwk.coloring import ccolors
+from fwk.coloring import ccolors
def main():
# timer
@@ -39,31 +37,38 @@ def main():
tms['total'].start()
# define flow variables
- alpha = 3*np.pi/180
- M_inf = 0.3
- c_ref = 1
- span = 1
+ alpha = 1*np.pi/180
+ M_inf = 0.72
dim = 3
+ # define dimension and mesh size
+ spn = 0.1 # wing span
+ lgt = 6.0 # channel length
+ hgt = 6.0 # channel height
+ wdt = 3.0 # channel width
+ c_ref = 1.0 # reference length
+ S_ref = spn # reference area
+ fms = 1.0 # farfield mesh size
+ nms = 0.01 # nearfield mesh size
+
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'spn' : span, 'lgt' : 3, 'wdt' : 3, 'hgt' : 3, 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
- # create mesh deformation handler
+ tms["msh"].start()
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
+ msh, gmshWriter = floD.mesh(dim, 'models/rae_25.geo', pars, ['field', 'wing', 'symmetry', 'downstream'])
morpher = floD.morpher(msh, dim, ['wing'])
- tms['msh'].stop()
+ tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing')
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
- solver.run()
+ cnvrgd = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
@@ -83,7 +88,6 @@ def main():
dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]))
dClX = np.sqrt(dClX)
dCdX = np.sqrt(dCdX)
-
# recover gradient wrt to AoA from mesh gradient
drot = np.array([[-np.sin(alpha), 0, np.cos(alpha)], [0, 0, 0], [-np.cos(alpha), 0, -np.sin(alpha)]])
dClAoA = 0
@@ -95,8 +99,8 @@ def main():
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' M alpha Cl Cd dCl_da dCd_da')
- print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alpha*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa))
+ print(' M alpha Cl Cd Cm')
+ print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/np.pi, solver.Cl, solver.Cd, solver.Cm))
# display timers
tms['total'].stop()
@@ -104,25 +108,25 @@ def main():
print('CPU statistics')
print(tms)
+ # visualize solution
+ floD.initViewer(pbl)
+
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ if (not cnvrgd):
+ raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
- if M_inf == 0.3 and alpha == 3*np.pi/180 and span == 1:
- tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 2.5, 5e-2)) # FD 2.526 (step 1e-5)
- tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.146, 5e-2)) # FD 0.1460 (step 1e-5)
- tests.add(CTest('dCl_dMsh', dClX, 2.016, 5e-2))
- tests.add(CTest('dCd_dMsh', dCdX, 0.133, 5e-2))
- tests.add(CTest('dCl/dAoA (msh)', dClAoA, 2.5, 5e-2))
- tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.146, 5e-2))
- elif M_inf == 0.7 and alpha == 5*np.pi/180 and span == 2:
- tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 4.9, 5e-2))
- tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.454, 5e-2))
- tests.add(CTest('dCl_dMsh', dClX, 2.571, 5e-2))
- tests.add(CTest('dCd_dMsh', dCdX, 0.204, 5e-2))
- tests.add(CTest('dCl/dAoA (msh)', dClAoA, 4.9, 5e-2))
- tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.454, 5e-2))
- else:
- raise Exception('Test not defined for this flow')
+ tests.add(CTest('iteration count', solver.nIt, 6, 1, forceabs=True))
+ tests.add(CTest('CL', solver.Cl, 0.61, 5e-2)) # 2D, 0.62
+ tests.add(CTest('CD', solver.Cd, 0.0016, 1e-2)) # 2D 0.0009 (3D refined 0.0006)
+ tests.add(CTest('CS', solver.Cs, 0.0000, 1e-3, forceabs=True))
+ tests.add(CTest('CM', solver.Cm, -0.116, 5e-2)) # 2D -0.117
+ tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 12.7, 5e-2)) # 2D 12.7, FD 12.65 (1e-6)
+ tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.073, 5e-2)) # 2D 0.051, FD, 0.073 (1e-6)
+ tests.add(CTest('dCl_dMsh', dClX, 38.0, 5e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 2.0, 5e-2))
+ tests.add(CTest('dCl/dAoA (msh)', dClAoA, 12.8, 5e-2)) # 2D 12.9
+ tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.073, 5e-2)) # 2D 0.051
tests.run()
# eof
diff --git a/dart/tests/rae_3.py b/dart/tests/rae_3.py
new file mode 100644
index 0000000..74a4120
--- /dev/null
+++ b/dart/tests/rae_3.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute the nonlinear flow around rectangular 3D rae2822
+# Adrien Crovato
+#
+# Test the direct and adjoint solvers, and the mesh morpher in 3D
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import numpy as np
+import dart.default as floD
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ alpha = 0
+ M_inf = 0.8
+ dim = 3
+
+ # define dimension and mesh size
+ spn = 1.0 # wing span
+ lgt = 6.0 # channel length
+ hgt = 6.0 # channel height
+ wdt = 3.0 # channel width
+ c_ref = 1.0 # reference length
+ S_ref = spn # reference area
+ fms = 1.0 # farfield mesh size
+ nms = 0.02 # nearfield mesh size
+
+ # parameters for mesh deformation
+ alfa = 1*np.pi/180
+ xc = 0.25
+ dz_max = 0.2*spn
+
+ # mesh the geometry and create mesh deformation handler
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms["msh"].start()
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
+ msh, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ morpher = floD.morpher(msh, dim, ['wing'])
+ tms["msh"].stop()
+
+ # set the problem and add fluid, initial/boundary conditions
+ tms['pre'].start()
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
+ solver = floD.newton(pbl) # create the solver now so element memory is up to date
+ tms['pre'].stop()
+
+ # deform the mesh (dummy rotation + bending)
+ print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET)
+ tms['deform'].start()
+ morpher.savePos()
+ rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]])
+ vNods = {}
+ for e in bnd.groups[0].tag.elems:
+ for n in e.nodes:
+ if n.no not in vNods:
+ vNods[n.no] = n
+ for no in vNods:
+ nod = np.matrix([[vNods[no].pos[0]-xc], [vNods[no].pos[2]]])
+ new = rMat * nod
+ vNods[no].pos[0] = new[0,0]
+ vNods[no].pos[2] = new[1,0] + dz_max/(spn*spn)*(vNods[no].pos[1]*vNods[no].pos[1])
+ # deform the mesh
+ morpher.deform()
+ gmshWriter.save(msh.name+"_def")
+ tms['deform'].stop()
+
+ # solve problem
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ cnvrgd = solver.run()
+ solver.save(gmshWriter)
+ tms['solver'].stop()
+
+ # solve adjoint problem
+ print(ccolors.ANSI_BLUE + 'PyGradient...' + ccolors.ANSI_RESET)
+ tms['adjoint'].start()
+ adjoint = floD.adjoint(solver, morpher)
+ adjoint.run()
+ adjoint.save(gmshWriter)
+ tms['adjoint'].stop()
+
+ # compute norm of gradient wrt to mesh
+ dClX = 0
+ dCdX = 0
+ for n in msh.nodes:
+ dClX += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]))
+ dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]))
+ dClX = np.sqrt(dClX)
+ dCdX = np.sqrt(dCdX)
+ # recover gradient wrt to AoA from mesh gradient
+ drot = np.array([[-np.sin(alfa), 0, np.cos(alfa)], [0, 0, 0], [-np.cos(alfa), 0, -np.sin(alfa)]])
+ dClAoA = 0
+ dCdAoA = 0
+ for n in bnd.nodes:
+ dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1], n.pos[2]]))
+ dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(dx)
+ dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(dx)
+
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' M alpha Cl Cd dCl_da dCd_da')
+ print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa))
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+
+ # visualize solution
+ floD.initViewer(pbl)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ if (not cnvrgd):
+ raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('iteration count', solver.nIt, 7, 1, forceabs=True))
+ tests.add(CTest('CL', solver.Cl, 0.21, 5e-2))
+ tests.add(CTest('CD', solver.Cd, 0.008, 1e-2))
+ tests.add(CTest('CS', solver.Cs, -0.020, 5e-2))
+ tests.add(CTest('CM', solver.Cm, -0.080, 1e-2))
+ tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 3.0, 5e-2)) # FD 3.0144 (1e-6)
+ tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.2, 5e-2)) # FD 0.2007 (1e-6)
+ tests.add(CTest('dCl_dMsh', dClX, 2.3, 5e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.2, 5e-2))
+ tests.add(CTest('dCl/dAoA (msh)', dClAoA, 3.0, 5e-2))
+ tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.2, 5e-2))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/validation/agard.py b/dart/validation/agard.py
index 3e048c0..86d5b9b 100644
--- a/dart/validation/agard.py
+++ b/dart/validation/agard.py
@@ -69,7 +69,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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 = 'wakeTip')
tms['pre'].stop()
# solve problem
diff --git a/dart/validation/onera.py b/dart/validation/onera.py
index 5739b25..6957094 100644
--- a/dart/validation/onera.py
+++ b/dart/validation/onera.py
@@ -67,7 +67,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
- pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, 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 = 'wakeTip')
tms['pre'].stop()
# solve problem
diff --git a/ext/amfe b/ext/amfe
index 87212cf..054ba00 160000
--- a/ext/amfe
+++ b/ext/amfe
@@ -1 +1 @@
-Subproject commit 87212cfb034a92cc83acfeea751f076a84c12779
+Subproject commit 054ba00c327e4f2745ff9c559832247f90c3cec9
--
GitLab