Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • Paul.Dechamps/dartflo
1 result
Show changes
Showing
with 1223 additions and 923 deletions
......@@ -26,7 +26,7 @@ def getParam():
# Arguments
args = parseargs()
p['Threads'] = args.k
p['Verb'] = args.verb
p['Verb'] = args.verb + 1
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/rae2822.geo') # Input file containing the model
p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model
......
......@@ -44,7 +44,7 @@ def getParam():
# Arguments
args = parseargs()
p['Threads'] = args.k
p['Verb'] = args.verb
p['Verb'] = args.verb + 1
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wbht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
......@@ -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)
......
......@@ -39,7 +39,7 @@ def getParam():
# Arguments
args = parseargs()
p['Threads'] = args.k
p['Verb'] = args.verb
p['Verb'] = args.verb + 1
# Input/Output
p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
......@@ -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)
......
......@@ -23,7 +23,7 @@ from fwk.wutils import parseargs
import dart
import tbox
def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
def mesh(dim, file, pars, bnds, wk = 'wake', wktp = None):
"""Initialize mesh and mesh writer
Parameters
......@@ -34,7 +34,7 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
file contaning mesh (w/o extention)
pars : dict
parameters for mesh
bnd : array of str
bnds : array of str
names of crack (wake) boundaries physical tag
wk : str (optional)
name of crack (wake) physical tag
......@@ -47,7 +47,8 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
# crack the mesh
mshCrck = tbox.MshCrack(msh, dim)
mshCrck.setCrack(wk)
mshCrck.addBoundaries(bnd)
for bnd in bnds:
mshCrck.addBoundary(bnd)
if dim == 3 and wktp is not None:
mshCrck.setExcluded(wktp)
mshCrck.run()
......@@ -57,7 +58,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 +98,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 +122,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)
......@@ -153,7 +151,7 @@ def picard(pbl):
problem formulation
"""
args = parseargs()
solver = dart.Picard(pbl, tbox.Gmres(), nthrds=args.k, vrb=args.verb+1)
solver = dart.Picard(pbl, tbox.Gmres(1, 1e-6, 30, 1e-8), nthrds=args.k, vrb=args.verb+2)
return solver
def newton(pbl):
......@@ -166,7 +164,7 @@ def newton(pbl):
"""
from tbox.solvers import LinearSolver
args = parseargs()
solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+1)
solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+2)
return solver
def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
......@@ -190,13 +188,14 @@ def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = '
name of physical tag contaning the symmetry plane (y)
"""
args = parseargs()
mshDef = tbox.MshDeform(msh, dim, nthrds=args.k)
mshDef = tbox.MshDeform(msh, tbox.Gmres(1, 1e-6, 30, 1e-8), dim, nthrds=args.k)
mshDef.setField(fld)
mshDef.addFixed(fxd)
for tag in fxd:
mshDef.addFixed(tag)
mshDef.addMoving(mov)
if dim == 3:
mshDef.setSymmetry(sym, 1)
mshDef.addInternal([wk, wk+'_'])
mshDef.addInternal(wk, wk+'_')
mshDef.initialize()
return mshDef
......@@ -210,8 +209,7 @@ def adjoint(solver, morpher):
solver : tbox.MshDeform object
mesh morpher
"""
from tbox.solvers import LinearSolver
adjoint = dart.Adjoint(solver, morpher, LinearSolver().pardiso())
adjoint = dart.Adjoint(solver, morpher)
return adjoint
def initViewer(problem):
......@@ -226,3 +224,20 @@ def initViewer(problem):
if not args.nogui:
from tbox.viewer import GUI
GUI().open( problem.msh.name )
def plot(x, y, cfg):
"""Plot data with check
Parameters
----------
x : array of float
x-data
y : array of float
y-data
cfg : dictionary
plot configuration
"""
args = parseargs()
if not args.nogui:
import tbox.utils as utils
utils.plot(x, y, cfg)
/* AGARD 445 wing */
// Initially generated by unsRgridWingGen.m
// For gmsh 4, use line 220 instead of line 221 (Bezier <- BSpline)
// Parameters
// domain and mesh
......@@ -344,7 +343,7 @@ MeshAlgorithm Surface{71,72,75,76} = 5;
/// 3D:
Mesh.Algorithm3D = 2;
Mesh.OptimizeNetgen = 1;
Mesh.Optimize = 1;
Mesh.Smoothing = 10;
Mesh.SmoothNormals = 1;
......@@ -354,7 +353,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};
......
......@@ -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
......
/*
* 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.
*/
/* NASA CRM - NASA Common Research Model
* CAD created by Guillem Battle I Capa
* Valid mesh obtained with gmsh 4.10.5
*/
/* Model information
* Length of fuselage: 62.79 m
* Chord length at wing root: 12.09 m
* Chord length at Yehudi break: 7.30 m
* Chord length at wing tip: 2.76 m
* Chord length at tail root: 5.45 m
* Chord length at tail tip: 2.33 m
* Reference wing area: 1/2 * 383.69 m^2
* Reference chord length: 7.01 m
*/
/// Constants
scl = 2; // scaling factor for lifting surfaces mesh size (testing purpose)
DefineConstant[ msFar = {36.8, Name "Mesh size on farfield" },
msFus = {0.628, Name "Mesh size on fuselage" },
msWingR = {0.121 * scl, Name "Mesh size on wing root" },
msWingY = {0.073 * scl, Name "Mesh size on wing Yehudi break" },
msWingT = {0.028 * scl, Name "Mesh size on wing tip" },
msTailR = {0.055 * scl, Name "Mesh size on tail root" },
msTailT = {0.023 * scl, Name "Mesh size on tail tip" } ];
/// Geometry definition
// Import
Geometry.OCCAutoFix = 1;
Geometry.OCCFixDegenerated = 1;
Geometry.OCCFixSmallEdges = 1;
Geometry.OCCFixSmallFaces = 1;
Geometry.OCCSewFaces = 1;
Geometry.Tolerance = 1e-5; // 0.01 mm
Geometry.OCCTargetUnit = "M";
Merge "crm.stp";
Coherence;
// Re-orient
ReverseMesh Surface{1:48}; //To reverse the direction of the mesh normal vectors
// Volume
Surface Loop(1) = {4:48};
Volume(1) = {1};
// Embedded
Line{3} In Surface{16};
Line{7} In Surface{16};
Line{15} In Surface{16};
Surface{1} In Volume{1};
Surface{2} In Volume{1};
Surface{3} In Volume{1};
// Physical groups
Physical Line("wakeExW") = {10,11,12,13};
Physical Line("teW") = {9,14};
Physical Line("wakeTipW") = {10};
Physical Line("wakeExT") = {1,2,5};
Physical Line("teT") = {6};
Physical Line("wakeTipT") = {1,2};
Physical Surface("symmetry") = {19,20,21,22};
Physical Surface("downstream") = {16};
Physical Surface("upstream") = {18};
Physical Surface("farfield") = {14,15,17};
Physical Surface("fuselage") = {13,23:45};
Physical Surface("wing") = {4,5,6,11,12};
Physical Surface("wing_") = {7,8,9,10};
Physical Surface("tail") = {47,48};
Physical Surface("tail_") = {46};
Physical Surface("wakeW") = {2,3};
Physical Surface("wakeT") = {1};
Physical Volume("field") = {1};
/// Mesh definition
// Farfield
Characteristic Length {34,35,36,37,38,39,40,43,44,45} = msFar; // farfield boundaries
Characteristic Length {3,4,7,8,15,41,42} = msFar; // wake
// Fuselage
Characteristic Length {11,12,13,51,61,50,73,74,62,63,64,65,47,29,70,46} = msFus; // body
Characteristic Length {51,52,53,54,55,56,59,60,89,85,86,90,91,92} = 0.5 * msFus; // nose and windshield
Characteristic Length {57,58,87,88} = 0.2 * msFus; // surface between windshield
Characteristic Length {49,72,28,33,32,24,31,30,75,71,48,76,69,47} = 0.75 * msFus; // wingbox
Characteristic Length {5,66,67,68} = 0.25 * msFus; // APU
Characteristic Length {46,11} = msFus; // Fuselage-wake intersection
// Wing
Characteristic Length {24, 25,26,27,14} = msWingR; // wing root
Characteristic Length {22, 20,21,23,9} = msWingY; // wing kink
Characteristic Length {18,10,16,19} = 3 * msWingT; // wing tip TE
Characteristic Length {17} = msWingT; // wing tip LE
// Tail
Characteristic Length {79:84,93:99, 100:102,77,78,103,6} = msTailR; // Tail root
Characteristic Length {1,2,105} = 2 * msTailT; // tail tip TE
Characteristic Length {104} = msTailT; // tail tip LE
// Mesh algos
Mesh.Algorithm = 6;
Mesh.Algorithm3D = 10;
Mesh.Optimize = 1;
Mesh.Smoothing = 10;
Mesh.SmoothNormals = 1;
Mesh 1;
Geometry.Tolerance = 1e-6;
Coherence Mesh;
File added
/* 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};
/* 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};
/* 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};
/* Onera M6 wing */
// Initially generated by unsRgridWingGen.m
// For gmsh 4, use line 425 instead of line 426 (Bezier <- BSpline)
// Parameters
// domain and mesh
......@@ -549,7 +548,7 @@ MeshAlgorithm Surface{71,72,75,76} = 5;
/// 3D:
Mesh.Algorithm3D = 2;
Mesh.OptimizeNetgen = 1;
Mesh.Optimize = 1;
Mesh.Smoothing = 10;
Mesh.SmoothNormals = 1;
......@@ -559,7 +558,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};
......
/* 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.Optimize = 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};
/* 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.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};
......@@ -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};
......
......@@ -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:
......
......@@ -69,12 +69,16 @@ class Adjoint;
class F0Ps;
class F0El;
class F0ElC;
class F0ElSpeedSound;
class F0ElRho;
class F0ElRhoL;
class F0ElRhoClamp;
class F0ElRhoLin;
class F0ElMach;
class F0ElMachL;
class F0ElMachClamp;
class F0ElMachLin;
class F0ElCp;
class F0ElCpL;
class F0ElCpClamp;
class F0ElCpLin;
class F1El;
class F1ElVi;
class F1Ct;
......@@ -86,4 +90,4 @@ class FaceEq;
} // namespace dart
#endif //DART_H
#endif // DART_H
......@@ -27,11 +27,14 @@
#include "wWakeResidual.h"
#include "wKuttaResidual.h"
#include "wLoadFunctional.h"
#include "wBlowing.h"
#include "wBlowingResidual.h"
#include "wMshData.h"
#include "wNode.h"
#include "wElement.h"
#include "wLinearSolver.h"
#include "wGmres.h"
#include "wMshDeform.h"
#include "wResults.h"
#include "wMshExport.h"
......@@ -51,12 +54,30 @@ using namespace dart;
/**
* @brief Initialize the adjoint solver
*/
Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol) : sol(_sol), morph(_morph), linsol(_linsol)
Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph) : sol(_sol), morph(_morph)
{
// Default parameters
nthreads = sol->nthreads;
verbose = sol->verbose;
// Linear solvers (if GMRES, use the same, but with a thighter tolerance)
if (sol->linsol->type() == LSolType::GMRES_ILUT)
{
std::shared_ptr<Gmres> gmres = std::make_shared<Gmres>(*dynamic_cast<Gmres *>(sol->linsol.get()));
gmres->setTolerance(1e-8);
alinsol = gmres;
}
else
alinsol = sol->linsol;
if (morph->linsol->type() == LSolType::GMRES_ILUT)
{
std::shared_ptr<Gmres> gmres = std::make_shared<Gmres>(*dynamic_cast<Gmres *>(morph->linsol.get()));
gmres->setTolerance(1e-12);
mlinsol = gmres;
}
else
mlinsol = morph->linsol;
// Update element gradients
sol->pbl->initGradElems();
......@@ -67,6 +88,13 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
dCdMsh.resize(sol->pbl->msh->nodes.size(), Eigen::Vector3d::Zero());
dClAoa = 0;
dCdAoa = 0;
// Display solver parameters
std::cout << "--- Adjoint solver ---\n"
<< "Aero inner solver: " << *alinsol
<< "Mesh inner solver: " << *mlinsol
<< "Number of threads: " << nthreads << "\n"
<< std::endl;
}
/**
......@@ -81,10 +109,13 @@ void Adjoint::run()
// Init
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Display solver parameters
std::cout << "--- Adjoint solver ---\n"
<< "Number of threads: " << nthreads << "\n"
<< std::endl;
// Display current freestream conditions
if (verbose > 0)
std::cout << std::fixed << std::setprecision(2)
<< "Computing gradients for Mach " << sol->pbl->M_inf << ", "
<< sol->pbl->alpha * 180 / 3.14159 << "deg AoA, "
<< sol->pbl->beta * 180 / 3.14159 << "deg AoS"
<< std::endl;
// Compute partial gradients of flow residuals and solve flow adjoint
tms["0-AdjFlo"].start();
......@@ -117,8 +148,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)
{
......@@ -165,7 +196,7 @@ void Adjoint::run()
<< std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lamCdMsh.data(), lamCdMsh.size()).norm() << std::endl;
// Display timers
if (verbose > 1)
if (verbose > 2)
std::cout << "Adjoint solver CPU" << std::endl
<< tms;
std::cout << std::endl;
......@@ -194,7 +225,7 @@ std::vector<double> Adjoint::solveMesh(std::vector<double> const &dX)
Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
Eigen::VectorXd dX_ = Eigen::Map<const Eigen::VectorXd>(dX.data(), dX.size());
// Solve
linsol->compute(dRx_X.transpose(), Eigen::Map<Eigen::VectorXd>(dX_.data(), dX_.size()), dR_);
mlinsol->compute(dRx_X.transpose(), Eigen::Map<Eigen::VectorXd>(dX_.data(), dX_.size()), dR_);
return std::vector<double>(dR.data(), dR.data() + dR.size());
}
......@@ -242,7 +273,7 @@ std::vector<double> Adjoint::solveFlow(std::vector<double> const &dU)
Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
Eigen::VectorXd dU_ = Eigen::Map<const Eigen::VectorXd>(dU.data(), dU.size());
// Solve
linsol->compute(dRu_U.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), dR_);
alinsol->compute(dRu_U.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), dR_);
return std::vector<double>(dR.data(), dR.data() + dR.size());
}
......@@ -332,18 +363,7 @@ std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsFlow()
// Compute gradients of lift and drag
Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size());
for (auto bnd : sol->pbl->bodies)
{
// gradients of loads
Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->msh->nodes.size());
this->buildGradientLoadsFlow(dL, *bnd);
// project to lift and drag directions
Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
for (size_t i = 0; i < bnd->nodes.size(); ++i)
{
dCl += dirL * dL.middleRows(i * 3, 3);
dCd += dirD * dL.middleRows(i * 3, 3);
}
}
this->buildGradientCoefficientsFlow(dCl, dCd, *bnd);
return std::vector<std::vector<double>>{std::vector<double>(dCl.data(), dCl.data() + dCl.size()), std::vector<double>(dCd.data(), dCd.data() + dCd.size())};
}
......@@ -369,23 +389,12 @@ std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsMesh()
// Compute the gradients of the lift and drag
Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
for (auto bnd : sol->pbl->bodies)
{
// gradients of loads
Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
this->buildGradientLoadsMesh(dL, *bnd);
// project to lift and drag directions
Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
for (size_t i = 0; i < bnd->nodes.size(); ++i)
{
dCl += dirL * dL.middleRows(i * 3, 3);
dCd += dirD * dL.middleRows(i * 3, 3);
}
}
this->buildGradientCoefficientsMesh(dCl, dCd, *bnd);
return std::vector<std::vector<double>>{std::vector<double>(dCl.data(), dCl.data() + dCl.size()), std::vector<double>(dCd.data(), dCd.data() + dCd.size())};
}
/**
* @brief Build the gradients of the loads with respect to the flow variables one a given body
* @brief Build gradients of the loads with respect to the flow variables on a given body
*/
void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd)
{
......@@ -404,7 +413,7 @@ void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
int rowi = bnd.nMap.at(e->nodes[i]);
int rowi = static_cast<int>(bnd.nMap.at(e->nodes[i]));
for (int m = 0; m < 3; ++m)
{
for (size_t j = 0; j < eV->nodes.size(); ++j)
......@@ -421,7 +430,43 @@ void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor
}
/**
* @brief Compute gradients of the residual with respect to angle of attack
* @brief Build gradients of the load coefficients with respect to the flow variables on a given body
*
* Note: this is essentially the same as buildGradientLoadsFlow, but with an additional projection performed inside the loop for performance
*/
void Adjoint::buildGradientCoefficientsFlow(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd)
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
// Build gradients
Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
tbb::parallel_for_each(bnd.groups[0]->tag->elems.begin(), bnd.groups[0]->tag->elems.end(), [&](Element *e) {
// Associated volume element
Element *eV = bnd.svMap.at(e);
// Build
Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->fluid->cP);
// Project
Eigen::RowVectorXd dCle = Eigen::RowVectorXd::Zero(eV->nodes.size()), dCde = Eigen::RowVectorXd::Zero(eV->nodes.size());
for (size_t i = 0; i < e->nodes.size(); ++i)
{
dCle += dirL * Ae.middleRows(i * 3, 3);
dCde += dirD * Ae.middleRows(i * 3, 3);
}
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
int rowj = eV->nodes[j]->row;
dCl(rowj) += dCle(j);
dCd(rowj) += dCde(j);
}
});
}
/**
* @brief Build gradients of the residual with respect to angle of attack
*/
void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
{
......@@ -451,7 +496,7 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
}
/**
* @brief Compute gradients of the lift and drag with respect to angle of attack
* @brief Build gradients of the lift and drag with respect to angle of attack
*/
void Adjoint::buildGradientLoadsAoa(double &dCl, double &dCd)
{
......@@ -466,7 +511,7 @@ void Adjoint::buildGradientLoadsAoa(double &dCl, double &dCd)
}
/**
* @brief Compute gradients of the residual with respect to mesh coordinates
* @brief Build gradients of the residual with respect to mesh coordinates
*/
void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
{
......@@ -482,7 +527,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
// Current element
Element *e = p.first;
//Upwind element
// Upwind element
Element *eU = p.second[0];
// Build elementary matrices
Eigen::MatrixXd Ae1, Ae2;
......@@ -518,6 +563,27 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
});
// Farfield B.C. are not taken into account because they have no influence on the solution
// Blowing boundary
for (auto blow : sol->pbl->bBCs)
{
tbb::parallel_for_each(blow->uE.begin(), blow->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
// Build elementary matrices
Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < p.first->nodes.size(); ++i)
{
size_t rowi = p.first->nodes[i]->row;
for (size_t j = 0; j < p.first->nodes.size(); ++j)
{
int rowj = p.first->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, A(i, sol->pbl->nDim * j + n)));
}
}
});
}
// Wake
for (auto wake : sol->pbl->wBCs)
{
......@@ -585,7 +651,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
}
/**
* @brief Compute gradients of the loads with respect to mesh coordinates, on a given body
* @brief Build gradients of the loads with respect to mesh coordinates on a given body
*/
void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd)
{
......@@ -605,7 +671,7 @@ void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
int rowi = bnd.nMap.at(e->nodes[i]);
int rowi = static_cast<int>(bnd.nMap.at(e->nodes[i]));
for (int m = 0; m < 3; ++m)
{
// surface
......@@ -618,7 +684,7 @@ void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor
// volume
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
size_t rowj = eV->nodes[j]->row;
int rowj = eV->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
T.push_back(Eigen::Triplet<double>(3 * rowi + m, sol->pbl->nDim * rowj + n, Aev(3 * i + m, sol->pbl->nDim * j + n)));
}
......@@ -630,32 +696,79 @@ void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor
dL.makeCompressed();
}
/**
* @brief Build gradients of the loads with respect to mesh coordinates on a given body
*
* Note: this is essentially the same as buildGradientLoadsMesh, but with an additional projection performed inside the loop for performance
*/
void Adjoint::buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd)
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
// Build gradients
Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
tbb::parallel_for_each(bnd.groups[0]->tag->elems.begin(), bnd.groups[0]->tag->elems.end(), [&](Element *e) {
// Associated volume element
Element *eV = bnd.svMap.at(e);
// Build
Eigen::MatrixXd Aes, Aev;
std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->fluid->cP, sol->pbl->nDim);
// Project
Eigen::RowVectorXd dCles = Eigen::RowVectorXd::Zero(sol->pbl->nDim * e->nodes.size()), dCdes = Eigen::RowVectorXd::Zero(sol->pbl->nDim * e->nodes.size());
Eigen::RowVectorXd dClev = Eigen::RowVectorXd::Zero(sol->pbl->nDim * eV->nodes.size()), dCdev = Eigen::RowVectorXd::Zero(sol->pbl->nDim * eV->nodes.size());
for (size_t i = 0; i < e->nodes.size(); ++i)
{
dCles += dirL * Aes.middleRows(i * 3, 3);
dCdes += dirD * Aes.middleRows(i * 3, 3);
dClev += dirL * Aev.middleRows(i * 3, 3);
dCdev += dirD * Aev.middleRows(i * 3, 3);
}
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
// surface
for (size_t j = 0; j < e->nodes.size(); ++j)
{
int rowj = e->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
{
dCl(sol->pbl->nDim * rowj + n) += dCles(sol->pbl->nDim * j + n);
dCd(sol->pbl->nDim * rowj + n) += dCdes(sol->pbl->nDim * j + n);
}
}
// volume
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
int rowj = eV->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
{
dCl(sol->pbl->nDim * rowj + n) += dClev(sol->pbl->nDim * j + n);
dCd(sol->pbl->nDim * rowj + n) += dCdev(sol->pbl->nDim * j + n);
}
}
});
}
/**
* @brief Write the results
*/
void Adjoint::save(MshExport *mshWriter, int n)
void Adjoint::save(MshExport *mshWriter, std::string const &suffix)
{
// Write files
std::cout << "Saving files... " << std::endl;
// setup results
std::cout << "Saving files "
<< std::setprecision(2)
<< "(Mach " << sol->pbl->M_inf << ", "
<< sol->pbl->alpha * 180 / 3.14159 << " deg AoA, "
<< sol->pbl->beta * 180 / 3.14159 << " deg AoS)"
<< std::endl;
// save results
Results results;
results.scalars_at_nodes["lambdaClPhi"] = &lamClFlo;
results.scalars_at_nodes["lambdaCdPhi"] = &lamCdFlo;
results.vectors_at_nodes["gradClMsh"] = &dClMsh;
results.vectors_at_nodes["gradCdMsh"] = &dCdMsh;
// save (whole mesh and bodies)
if (n > 0)
{
mshWriter->save(sol->pbl->msh->name + "_adjoint_" + std::to_string(n), results);
for (auto bnd : sol->pbl->bodies)
bnd->save(bnd->groups[0]->tag->name + "adjoint_" + std::to_string(n), results);
}
else
{
mshWriter->save(sol->pbl->msh->name + "_adjoint", results);
for (auto bnd : sol->pbl->bodies)
bnd->save(bnd->groups[0]->tag->name + "_adjoint", results);
}
mshWriter->save(sol->pbl->msh->name + "_adjoint" + suffix, results);
}
void Adjoint::write(std::ostream &out) const
......
......@@ -37,9 +37,10 @@ namespace dart
class DART_API Adjoint : public fwk::wSharedObject
{
public:
std::shared_ptr<Newton> sol; ///< Newton solver
std::shared_ptr<tbox::MshDeform> morph; ///< mesh morpher
std::shared_ptr<tbox::LinearSolver> linsol; ///< linear solver
std::shared_ptr<Newton> sol; ///< Newton solver
std::shared_ptr<tbox::MshDeform> morph; ///< mesh morpher
std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
std::shared_ptr<tbox::LinearSolver> mlinsol; ///< linear solver (adjoint mesh)
int nthreads; ///< number of threads for the assembly
int verbose; ///< display more info
......@@ -61,11 +62,11 @@ private:
fwk::Timers tms; ///< internal timers
public:
Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol);
Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph);
virtual ~Adjoint() { std::cout << "~Adjoint()\n"; }
void run();
void save(tbox::MshExport *mshWriter, int n = 0);
void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
// Partial gradients of the mesh residuals
void linearizeMesh();
......@@ -87,19 +88,27 @@ public:
std::vector<std::vector<double>> computeGradientCoefficientsMesh();
#ifndef SWIG
// Getters for partials (required by VII-BLASTER)
Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_U() const { return dRu_U; }
Eigen::VectorXd getRu_A() const { return dRu_A; }
Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_X() const { return dRu_X; }
Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
virtual void write(std::ostream &out) const override;
#endif
private:
// Partial gradients wrt with respect to flow variables
void buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd);
void buildGradientCoefficientsFlow(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd);
// Partial gradients with respect to angle of attack
void buildGradientResidualAoa(Eigen::VectorXd &dR);
void buildGradientLoadsAoa(double &dCl, double &dCd);
// Partial gradients with respect to mesh coordinates
void buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
void buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd);
void buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd);
};
} // namespace dart
#endif //WADJOINT_H
#endif // WADJOINT_H
......@@ -77,9 +77,16 @@ void Initial::write(std::ostream &out) const
Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta, bool pin) : Assign(_msh, no)
{
f = new F0PsPhiInf(dim, alpha, beta);
// Only retain the first node, so that the DOF associated to this node will be pinned
// Only retain one node, so that the DOF associated to this node will be pinned
if (pin)
this->nodes.resize(1);
{
// Pin node with min x min y min z
auto it = std::min_element(nodes.begin(), nodes.end(), [](Node *n1, Node *n2) {
return n1->pos[0] < n2->pos[0] || (n1->pos[0] == n2->pos[0] && n1->pos[1] < n2->pos[1]) || (n1->pos[0] == n2->pos[0] && n1->pos[1] == n2->pos[1] && n1->pos[2] < n2->pos[2]);
});
nodes.resize(1);
nodes[0] = *it;
}
}
Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name)
{
......