Commit 8c6c39d2 authored by 's avatar

Merge remote-tracking branch 'upstream/master' into fix/inverted-triangles

parents bc34ad0d e38421d7
Changelog
=========
v4.5.0
v5.0.0
------
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
- Add periodic boundary condition meshing in 2D (!733)
- Adjust boundary layer thickness in corners in 2D (!739)
- Add non-O BL meshing in 2D (!757)
- Fix issue with reading CCM files due to definition of default arrays rather than a vector (!797)
**Library**
- Added in sum factorisation version for pyramid expansions and orthogonal
......@@ -152,6 +155,7 @@ v4.4.0
(!712)
- 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720)
- Add basic gmsh cad (.geo) reader to the meshing system (!731)
- Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boundary layers in 2D (!736)
- Change variable names in mcf file to make more sense (!736)
......
......@@ -308,10 +308,6 @@ IF( NEKTAR_USE_TINYXML_STL )
ADD_DEFINITIONS( -DTIXML_USE_STL)
ENDIF( NEKTAR_USE_TINYXML_STL )
IF( NEKTAR_USE_DIRECT_BLAS_CALLS )
ADD_DEFINITIONS(-DNEKTAR_USING_DIRECT_BLAS_CALLS)
ENDIF( NEKTAR_USE_DIRECT_BLAS_CALLS )
IF( NEKTAR_USE_EXPRESSION_TEMPLATES )
ADD_DEFINITIONS(-DNEKTAR_USE_EXPRESSION_TEMPLATES -DNEKTAR_USING_CMAKE)
ENDIF( NEKTAR_USE_EXPRESSION_TEMPLATES )
......
......@@ -437,6 +437,7 @@ Compressible flow is characterised by abrupt changes in density within the flow
\begin{equation}\label{eq:sensor}
S_e=\frac{||\rho^p_e-\rho^{p-1}_e||_{L_2}}{||\rho_e^p||_{L_2}}
\end{equation}
by default the comparison is made with the $p-1$ solution, but this can be changed by setting the parameter \inltt{SensorOffset}.
An artificial diffusion term is introduced locally to the Euler equations to deal with flow discontinuity and the consequential numerical oscillations. Two models are implemented, a non-smooth and a smooth artificial viscosity model.
\subsubsection{Non-smooth artificial viscosity model}
For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added laplacian term on right hand side of equation \ref{eq:euler}. The diffusivity of the system is controlled by a variable viscosity coefficient $\epsilon$. The value of $\epsilon$ is dependent on $\epsilon_0$, which is the maximum viscosity that is dependent on the polynomial order ($p$), the mesh size ($h$) and the maximum wave speed and the local sensor value. Based on pre-defined sensor threshold values, the variable viscosity is set accordingly
......@@ -512,4 +513,53 @@ The polynomial order in each element can be adjusted based on the sensor value t
\end{equation}
For now, the threshold values $s_e$, $s_{ds}$, $s_{sm}$ and $s_{fl}$ are determined empirically by looking at the sensor distribution in the domain. Once these values are set, two .txt files are outputted, one that has the composites called VariablePComposites.txt and one with the expansions called VariablePExpansions.txt. These values have to copied into a new .xml file to create the adapted mesh.
\subsection{Quasi-1D nozzle flow}
A quasi-1D inviscid flow (flow with area variation) can be obtained using the
\inltt{Quasi1D} forcing in a 1D solution of the Euler equations:
\begin{lstlisting}[style=XMLStyle]
<FORCING>
<FORCE TYPE="Quasi1D">
<AREAFCN> Area </AREAFCN>
</FORCE>
</FORCING>
\end{lstlisting}
in this case a function named \inltt{Area} must be specified in the \inltt{CONDITIONS} section of the session file.
In this case, it is possible to prescribe the inflow conditions in terms of stagnation properties (density and pressure)
by using the following boundary condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="rho" USERDEFINEDTYPE="StagnationInflow" VALUE="rhoStag" />
<D VAR="rhou" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="E" USERDEFINEDTYPE="StagnationInflow" VALUE="pStag/(Gamma-1)" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
\subsection{Axi-symmetric flow}
An axi-symmetric inviscid flow (with symmetry axis on x=0) can be obtained using
the \inltt{AxiSymmetric} forcing in a 2D solution of the Euler equations:
\begin{lstlisting}[style=XMLStyle]
<FORCING>
<FORCE TYPE="AxiSymmetric">
</FORCE>
</FORCING>
\end{lstlisting}
The \inltt{StagnationInflow} boundary condition can also be used in this case.
Also, by defining the geometry with \inltt{<GEOMETRY DIM="2" SPACE="3">} (i.e. a two-dimensional
mesh in three-dimensional space) and adding the \inltt{rhow} variable, we obtain an axi-symmetric
flow with swirl, in which case the \inltt{StagnationInflow} boundary condition allows prescribing \inltt{rhow}:
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="rho" USERDEFINEDTYPE="StagnationInflow" VALUE="rhoStag" />
<D VAR="rhou" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="rhov" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="rhow" USERDEFINEDTYPE="StagnationInflow" VALUE="x" />
<D VAR="E" USERDEFINEDTYPE="StagnationInflow" VALUE="pStag/(Gamma-1)" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
This diff is collapsed.
......@@ -842,22 +842,22 @@ example which generates a 2D NACA wing.
<MESHING>
<INFORMATION>
<I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2D" />
<I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2D" />
</INFORMATION>
<PARAMETERS>
<P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" />
<P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" />
<!-- 2D Domain !-->
<P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" />
<P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" />
</PARAMETERS>
</MESHING>
......@@ -867,10 +867,11 @@ In all cases the mesh generator needs two pieces of information and four
parameters. It firstly needs to know the CAD file with which to work. In the
example above this is listed as a 4 digit number, this is because the mesh
generator is equiped with a NACA wing generator. In all other cases this
parameter would be a STEP file. Secondly, what type of mesh to make, the options
are \inltt{EULER} and \inltt{BL} for 3D meshes and \inltt{2D} and \inltt{2DBL}
parameter would be the name of a CAD file (in either STEP or GEO format).
Secondly, what type of mesh to make, the options
are \inltt{EULER} and \inltt{BndLayer} for 3D meshes and \inltt{2D} and \inltt{2DBndLayer}
for 2D meshes. In the case of \inltt{EULER} the mesh will be made with only
tetrahedra. For \inltt{BL} the mesh generator will attempt to insert a single
tetrahedra. For \inltt{BndLayer} the mesh generator will attempt to insert a single
macro prism layer onto the geometry surface. This option requires additional
parameters. This is similar for the 2D scenarios. The automatic mesh
specification system requires three parameters to build the specification of a
......@@ -889,35 +890,35 @@ An example is shown.
<MESHING>
<INFORMATION>
<I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2DBL" />
<I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2DBndLayer" />
</INFORMATION>
<PARAMETERS>
<P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" />
<P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" />
<!-- Boundary layer !-->
<P PARAM="BLSurfs" VALUE="5,6" />
<P PARAM="BLThick" VALUE="0.03" />
<P PARAM="BLLayers" VALUE="4" />
<P PARAM="BLProg" VALUE="2.0" />
<P PARAM="BndLayerSurfaces" VALUE="5,6" />
<P PARAM="BndLayerThickness" VALUE="0.03" />
<P PARAM="BndLayerLayers" VALUE="4" />
<P PARAM="BndLayerProgression" VALUE="2.0" />
<!-- 2D Domain !-->
<P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" />
<P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" />
</PARAMETERS>
</MESHING>
</NEKTAR>
\end{lstlisting}
A list of the CAD surfaces which will have a prism generated on is described by
\inltt{BLSurfs} and the thickness of the boundary to aim for is \inltt{BLThick}.
\inltt{BndLayerSurfaces} and the thickness of the boundary to aim for is \inltt{BndLayerThickness}.
%
The mesh generator has been created with a range of error messages to aid in
debugging. If you encounter an error and the mesh generator fails, run \nm with
......@@ -925,6 +926,32 @@ the verbose \inltt{-v} flag and send the stdout with the .mcf and .step files
to \inltt{m.turner14@imperial.ac.uk}. Without the feedback this functionality
cannot improve.
\subsubsection{GEO format}
Recent developments have been made to facilitate the generation of meshes from
simple 2D geometries. The GEO file format, used by Gmsh, is a popular option
that allows the user to script geometrical and meshing operations without the
need of a GUI. A simplified reader has been implemented in NekMesh for 2D geometries.
Although very basic this reader may be extended in the future to cover a wider
range of geometrical features.
For a full description of the GEO format the user should refer to Gmsh's
documentation. The following commands are currently supported:
\begin{itemize}
\item \inltt{//} (comments)
\item \inltt{Point}
\item \inltt{Line}
\item \inltt{Spline}
\item \inltt{Line Loop}
\item \inltt{Plane Surface}
\end{itemize}
At the present time, NekMesh does not support the full scripting capabilities of the
GEO format. The used GEO files should be a straightforward succession of entity
creations (see list above). This should however allow for the creation of quite
a wide range of 2D geometries by transformation of arbitrary curves into generic
splines.
%%% Local Variables:
%%% mode: latex
......
......@@ -69,9 +69,7 @@ As an example, consider:
\end{lstlisting}
This will create a sequence of files named \inltt{MyFile\_*\_fc.vtu} containing isocontours.
The result will be output every 100 time steps. Output directly to
\inltt{.vtu} or \inltt{.dat} is currently only supported for isocontours.
In other cases, the output should be a \inltt{.fld} file.
The result will be output every 100 time steps.
\subsection{Time-averaged fields}
......
......@@ -147,7 +147,7 @@ class IProductWRTDerivBase_StdMat : public Operator
{
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim();
m_coordim = pCollExp[0]->GetCoordim();
int nqtot = m_stdExp->GetTotPoints();
int nmodes = m_stdExp->GetNcoeffs();
......
......@@ -8,6 +8,7 @@ SET(FieldUtilsHeaders
InputModules/InputPts.h
InputModules/InputNek5000.h
InputModules/InputSemtex.h
OutputModules/OutputFileBase.h
OutputModules/OutputInfo.h
OutputModules/OutputTecplot.h
OutputModules/OutputVtk.h
......@@ -20,6 +21,7 @@ SET(FieldUtilsHeaders
ProcessModules/ProcessBoundaryExtract.h
ProcessModules/ProcessCombineAvg.h
ProcessModules/ProcessConcatenateFld.h
ProcessModules/ProcessCreateExp.h
ProcessModules/ProcessDeform.h
ProcessModules/ProcessDisplacement.h
ProcessModules/ProcessEquiSpacedOutput.h
......@@ -30,6 +32,7 @@ SET(FieldUtilsHeaders
ProcessModules/ProcessInterpField.h
ProcessModules/ProcessInterpPoints.h
ProcessModules/ProcessInterpPointDataToFld.h
ProcessModules/ProcessInterpPtsToPts.h
ProcessModules/ProcessIsoContour.h
ProcessModules/ProcessJacobianEnergy.h
ProcessModules/ProcessMapping.h
......@@ -57,6 +60,7 @@ SET(FieldUtilsSources
InputModules/InputPts.cpp
InputModules/InputNek5000.cpp
InputModules/InputSemtex.cpp
OutputModules/OutputFileBase.cpp
OutputModules/OutputInfo.cpp
OutputModules/OutputTecplot.cpp
OutputModules/OutputVtk.cpp
......@@ -69,6 +73,7 @@ SET(FieldUtilsSources
ProcessModules/ProcessBoundaryExtract.cpp
ProcessModules/ProcessCombineAvg.cpp
ProcessModules/ProcessConcatenateFld.cpp
ProcessModules/ProcessCreateExp.cpp
ProcessModules/ProcessDeform.cpp
ProcessModules/ProcessDisplacement.cpp
ProcessModules/ProcessEquiSpacedOutput.cpp
......@@ -79,6 +84,7 @@ SET(FieldUtilsSources
ProcessModules/ProcessInterpField.cpp
ProcessModules/ProcessInterpPoints.cpp
ProcessModules/ProcessInterpPointDataToFld.cpp
ProcessModules/ProcessInterpPtsToPts.cpp
ProcessModules/ProcessIsoContour.cpp
ProcessModules/ProcessJacobianEnergy.cpp
ProcessModules/ProcessMapping.cpp
......
......@@ -50,19 +50,10 @@
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/DisContField1D.h>
#include <MultiRegions/DisContField2D.h>
#include <MultiRegions/DisContField3D.h>
#include <MultiRegions/DisContField3DHomogeneous1D.h>
#include <MultiRegions/DisContField3DHomogeneous2D.h>
#include "FieldUtilsDeclspec.h"
using namespace std;
namespace Nektar
{
namespace FieldUtils
......@@ -74,8 +65,7 @@ struct Field
: m_verbose(false), m_declareExpansionAsContField(false),
m_declareExpansionAsDisContField(false),
m_requireBoundaryExpansion(false), m_writeBndFld(false),
m_fldToBnd(false), m_addNormals(false),
m_setUpEquiSpacedFields(false), m_fieldPts(LibUtilities::NullPtsField)
m_addNormals(false), m_fieldPts(LibUtilities::NullPtsField)
{
}
......@@ -87,9 +77,12 @@ struct Field
}
}
bool m_verbose;
vector<LibUtilities::FieldDefinitionsSharedPtr> m_fielddef;
vector<vector<double> > m_data;
vector<MultiRegions::ExpListSharedPtr> m_exp;
std::vector<LibUtilities::FieldDefinitionsSharedPtr> m_fielddef;
std::vector<std::vector<double> > m_data;
std::vector<MultiRegions::ExpListSharedPtr> m_exp;
std::vector<std::string> m_variables;
int m_numHomogeneousDir;
bool m_declareExpansionAsContField;
bool m_declareExpansionAsDisContField;
......@@ -101,20 +94,14 @@ struct Field
LibUtilities::CommSharedPtr m_comm;
LibUtilities::SessionReaderSharedPtr m_session;
SpatialDomains::MeshGraphSharedPtr m_graph;
LibUtilities::PtsIOSharedPtr m_ptsIO;
map<string, vector<string> > m_inputfiles;
std::map<std::string, std::vector<std::string> > m_inputfiles;
bool m_writeBndFld;
vector<unsigned int> m_bndRegionsToWrite;
bool m_fldToBnd;
std::vector<unsigned int> m_bndRegionsToWrite;
bool m_addNormals;
bool m_setUpEquiSpacedFields;
LibUtilities::PtsFieldSharedPtr m_fieldPts;
MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap;
FIELD_UTILS_EXPORT MultiRegions::ExpListSharedPtr SetUpFirstExpList(
......@@ -418,10 +405,9 @@ struct Field
* @return Reader for @p filename.
*/
FIELD_UTILS_EXPORT LibUtilities::FieldIOSharedPtr FieldIOForFile(
string filename)
std::string filename)
{
LibUtilities::CommSharedPtr c = m_session ? m_session->GetComm() :
LibUtilities::GetCommFactory().CreateInstance("Serial", 0, 0);
LibUtilities::CommSharedPtr c = m_comm;
string fmt = LibUtilities::FieldIO::GetFileType(filename, c);
map<string, LibUtilities::FieldIOSharedPtr>::iterator it =
m_fld.find(fmt);
......@@ -440,7 +426,9 @@ struct Field
}
FIELD_UTILS_EXPORT MultiRegions::ExpListSharedPtr AppendExpList(
int NumHomogeneousDir, string var = "DefaultVar", bool NewField = false)
int NumHomogeneousDir,
std::string var = "DefaultVar",
bool NewField = false)
{
if (var.compare("DefaultVar") == 0 && m_requireBoundaryExpansion)
{
......@@ -691,8 +679,6 @@ struct Field
tmp = MemoryManager<MultiRegions::ContField3D>::
AllocateSharedPtr(*tmp2, m_graph, var);
m_locToGlobalMap = tmp2->GetLocalToGlobalMap();
}
}
else if (m_declareExpansionAsDisContField)
......@@ -731,10 +717,21 @@ struct Field
return tmp;
}
FIELD_UTILS_EXPORT void ClearField()
{
m_session = LibUtilities::SessionReaderSharedPtr();
m_graph = SpatialDomains::MeshGraphSharedPtr();
m_fieldPts = LibUtilities::NullPtsField;
m_exp.clear();
m_fielddef = std::vector<LibUtilities::FieldDefinitionsSharedPtr>();
m_data = std::vector<std::vector<NekDouble> > ();
m_variables.clear();
}
private:
/// Map to store FieldIO instances. Key is the reader type, value is the
/// FieldIO object.
map<string, LibUtilities::FieldIOSharedPtr> m_fld;
std::map<std::string, LibUtilities::FieldIOSharedPtr> m_fld;
};
typedef boost::shared_ptr<Field> FieldSharedPtr;
......
......@@ -53,7 +53,7 @@ ModuleKey InputDat::m_className[1] = {
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eInputModule, "dat"),
InputDat::create,
"Reads Tecplot dat file for FE block triangular format."),
"Reads Tecplot dat file for FE block triangular format.")
};
/**
......@@ -77,18 +77,8 @@ InputDat::~InputDat()
*/
void InputDat::Process(po::variables_map &vm)
{
if (m_f->m_verbose)
{
if (m_f->m_comm->TreatAsRankZero())
{
cout << "Processing input dat file" << endl;
}
}
string line, word, tag;
string line;
std::ifstream datFile;
stringstream s;
// Open the file stream.
string fname = m_f->m_inputfiles["dat"][0];
......@@ -101,7 +91,7 @@ void InputDat::Process(po::variables_map &vm)
}
// read variables
// currently assum there are x y and z coordinates
// currently assume there are x y and z coordinates
int dim = 3;
vector<string> fieldNames;
while (!datFile.eof())
......@@ -153,6 +143,9 @@ void InputDat::Process(po::variables_map &vm)
dim, fieldNames, pts);
m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
m_f->m_fieldPts->SetConnectivity(ptsConn);
// save field names
m_f->m_variables = fieldNames;
}
/**
......
......@@ -64,6 +64,16 @@ public:
return "InputDat";
}
virtual std::string GetModuleDescription()
{
return "Processing input dat file";
}
virtual ModulePriority GetModulePriority()
{
return eCreatePts;
}
private:
void ReadTecplotFEBlockZone(std::ifstream &datFile,
string &line,
......
......@@ -40,9 +40,6 @@ using namespace std;
#include "InputFld.h"
using namespace Nektar;
static std::string npts = LibUtilities::SessionReader::RegisterCmdLineArgument(
"NumberOfPoints", "n", "Define number of points to dump output");
namespace Nektar
{
namespace FieldUtils
......@@ -50,17 +47,18 @@ namespace FieldUtils
ModuleKey InputFld::m_className[4] = {
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eInputModule, "fld"), InputFld::create, "Reads Fld file."),
GetModuleFactory().RegisterCreatorFunction(ModuleKey(eInputModule, "chk"),
InputFld::create,
"Reads checkpoint file."),
GetModuleFactory().RegisterCreatorFunction(ModuleKey(eInputModule, "rst"),
InputFld::create,
"Reads restart file."),
ModuleKey(eInputModule, "fld"),
InputFld::create, "Reads Fld file."),
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eInputModule, "chk"),
InputFld::create, "Reads checkpoint file."),
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eInputModule, "rst"),
InputFld::create, "Reads restart file."),
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eInputModule, "bse"),
InputFld::create,
"Reads stability base-flow file.")};
InputFld::create, "Reads stability base-flow file.")
};
/**
* @brief Set up InputFld object.
......@@ -86,170 +84,63 @@ InputFld::~InputFld()
*/
void InputFld::Process(po::variables_map &vm)
{
int i;
string fileName = m_config["infile"].as<string>();
if (m_f->m_verbose)
{
if (m_f->m_comm->TreatAsRankZero())
{
cout << "Processing input fld file" << endl;
}
}
int i, j;
string fldending;
// Determine appropriate field input
if (m_f->m_inputfiles.count("fld") != 0)
{
fldending = "fld";
}
else if (m_f->m_inputfiles.count("chk") != 0)
{
fldending = "chk";
}
else if (m_f->m_inputfiles.count("rst") != 0)
{
fldending = "rst";