Commit d63376fb authored by Kilian Lackhove's avatar Kilian Lackhove

Merge branch 'master' into fix/FieldConvert-list

parents 96954c5c d07c83a4
......@@ -24,10 +24,12 @@ v5.0.0
- Fix issue with reading CCM files due to definition of default arrays
rather than a vector (!797)
- Fix inverted triangles and small memory issue in surface meshing (!798)
- Additional curve types in GEO reader: BSpline, Circle, Ellipse (!800)
**FieldConvert**:
- Add input module for Semtex field files (!777)
- Fixed interppoints module (!760)
- Move StreamFunction utility to a FieldConvert module (!809)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
......@@ -40,6 +42,9 @@ v4.4.2
**Library**
- Fix ability to set default implementation in Collections and added an option
to set eNoCollections in FieldConvert as default (!789)
- Fix performance issue in ProcessIsoContour in relation to memory consumption
(!821)
- Fix performance issue with ExtractPhysToBndElmt (!796)
**Utilities**
- Fix vtkToFld missing dependency which prevented compiling with VTK 7.1 (!808)
......
......@@ -60,8 +60,9 @@ IF (NEKTAR_USE_FFTW)
# Test if FFTW path is a system path. Only add to include path if not an
# implicitly defined CXX include path (due to GCC 6.x now providing its own
# version of some C header files and -isystem reorders include paths).
GET_FILENAME_COMPONENT(X ${CMAKE_CXX_IMPLICIT_INCLUDE_DIRECTORIES} ABSOLUTE)
GET_FILENAME_COMPONENT(X "${CMAKE_CXX_IMPLICIT_INCLUDE_DIRECTORIES}" ABSOLUTE)
GET_FILENAME_COMPONENT(Y ${FFTW_INCLUDE_DIR} ABSOLUTE)
IF (NOT Y MATCHES ".*${X}.*")
INCLUDE_DIRECTORIES(SYSTEM ${FFTW_INCLUDE_DIR})
ENDIF()
......
......@@ -475,3 +475,34 @@ volume = {163},
year = {2016}
}
@article{rospjo16,
title={Eigensolution analysis of spectral/hp continuous Galerkin approximations to advection--diffusion problems: Insights into spectral vanishing viscosity},
author={Moura, RC and Sherwin, SJ and Peir{\'o}, Joaquim},
journal={Journal of Computational Physics},
volume={307},
pages={401--422},
year={2016},
publisher={Elsevier}
}
@article{yvsiouei93,
title={Legendre pseudospectral viscosity method for nonlinear conservation laws},
author={Maday, Yvon and Kaber, Sidi M Ould and Tadmor, Eitan},
journal={SIAM Journal on Numerical Analysis},
volume={30},
number={2},
pages={321--342},
year={1993},
publisher={SIAM}
}
@article{rosh06,
title={Stabilisation of spectral/hp element methods through spectral vanishing viscosity: Application to fluid mechanics modelling},
author={Kirby, Robert M and Sherwin, Spencer J},
journal={Computer methods in applied mechanics and engineering},
volume={195},
number={23},
pages={3128--3144},
year={2006},
publisher={Elsevier}
}
\ No newline at end of file
......@@ -830,7 +830,6 @@ the advection term using the pressure inverse mass matrix. It can be used just i
<I PROPERTY="SmoothAdvection" VALUE="True"/>
\end{lstlisting}
\item \inltt{SpectralVanishingViscosity}: activates a stabilization technique
which increases the viscosity on the modes with the highest frequencies.
\begin{lstlisting}[style=XMLStyle]
......@@ -839,7 +838,34 @@ which increases the viscosity on the modes with the highest frequencies.
In a Quasi-3D simulation, this will affect both the Fourier and the spectral/hp expansions.
To activate them independently, use \inltt{SpectralVanishingViscositySpectralHP}
and \inltt{SpectralVanishingViscosityHomo1D}.
and \inltt{SpectralVanishingViscosityHomo1D}. \\
There are three spectral vanishing viscosity kernels available:
\begin{center}
\footnotesize
\begin{tabular}{lcc}
\toprule
{SVV Kernel} & {\texttt{SpectralVanishingViscosity}} \\
\midrule
Exponential Kernel & \texttt{True} \\
Power Kernel & \texttt{PowerKernel} \\
DG Kernel & \texttt{DGKernel} \\
\bottomrule
\end{tabular}
\end{center}
The Exponential kernel is based on the work of Maday et al. \cite{yvsiouei93},
its extension to 2D can be found in \cite{rosh06}. A diffusion coefficient can
be specified which defines the base magnitude of the viscosity; this parameter
is scaled by $h/p$. SVV viscosity is activated for expansion modes greater than
the product of the cut-off ratio and the expansion order. The Power kernel is a
smooth function with no cut-off frequency; it focusses on a narrower band of
higher expansion modes as the polynomial order increases. The cut-off ratio
parameter for the Power kernel corresponds to the power ratio, see Moura et al.
\cite{rospjo16}. The DG-Kernel is an attempt to match the dissipation of CG-SVV
to DG schemes of lower expansion orders. This kernel does not require any parameters
although the diffusion coefficient can still be modified.
\item \inltt{DEALIASING}: activates the 3/2 padding rule on the advection term
of a Quasi-3D simulation.
......@@ -856,7 +882,6 @@ stabilize the simulation. This method is based on the work of Kirby and Sherwin
\end{itemize}
\subsection{Parameters}
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file:
......@@ -869,7 +894,7 @@ the session file:
\item \inltt{MinSubSteps}: perform a minimum number of substeps in sub-stepping algorithm (default is 1)
\item \inltt{MaxSubSteps}: perform a maxmimum number of substeps in sub-stepping algorithm otherwise exit (default is 100)
\item \inltt{SVVCutoffRatio}: sets the ratio of Fourier frequency not affected by the SVV technique (default value = 0.75, i.e. the first 75\% of frequency are not damped)
\item \inltt{SVVDiffCoeff}: sets the SVV diffusion coefficient (default value = 0.1)
\item \inltt{SVVDiffCoeff}: sets the SVV diffusion coefficient (default value = 0.1 (Exponential and Power kernel), 1 (DG-Kernel))
\end{itemize}
\subsection{Womersley Boundary Condition}
......
......@@ -214,6 +214,7 @@ openany, % A chapter may start on either a recto or verso page.
}
\lstset{%
escapeinside={(*}{*)},%
breaklines=true,
}
\usepackage{tikz}
......
......@@ -178,6 +178,7 @@ possibly also Reynolds stresses) into single file;
\item \inltt{scalargrad}: Computes scalar gradient field;
\item \inltt{scaleinputfld}: Rescale input field by a constant factor;
\item \inltt{shear}: Computes time-averaged shear stress metrics: TAWSS, OSI, transWSS, TAAFI, TACFI, WSSG;
\item \inltt{streamfunction}: Calculates stream function of a 2D incompressible flow.
\item \inltt{surfdistance}: Computes height of a prismatic boundary layer mesh and projects onto the surface (for e.g. $y^+$ calculation).
\item \inltt{vorticity}: Computes the vorticity field.
\item \inltt{wss}: Computes wall shear stress field.
......@@ -543,7 +544,7 @@ In order to interpolate 1D data to a $n$D field, specify the matching coordinate
the output field using the \inltt{interpcoord} argument:
%
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m interppointdatatofld:frompts=1D-file1.pts:interppointdatatofld=1 \
FieldConvert -m interppointdatatofld:frompts=1D-file1.pts:interpcoord=1 \
3D-file1.xml 3D-file1.fld
\end{lstlisting}
%
......@@ -873,6 +874,26 @@ The argument \inltt{N} and \inltt{fromfld} are compulsory arguments that respect
The input \inltt{.fld} files are the outputs of the \textit{wss} module. If they do not contain the surface normals (an optional output of the \textit{wss} modle), then the \textit{shear} module will not compute the last metric, |WSSG|.
%
%
%
\subsection{Stream function of a 2D incompressible flow: \textit{streamfunction} module}
The streamfunction module calculates the stream function of a 2D incompressible flow, by
solving the Poisson equation
\[
\nabla^2 \psi = -\omega
\]
where $\omega$ is the vorticity. Note that this module applies the same boundary conditions
specified for the y-direction velocity component \inltt{v} to the stream function,
what may not be the most appropriate choice.
To use this module, the user can run
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m streamfunction test.xml test.fld test-streamfunc.fld
\end{lstlisting}
where the file \inltt{test-streamfunc.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}.
%
%
......
......@@ -938,10 +938,20 @@ 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{//} (i.e. comments)
\item \inltt{Point}
\item \inltt{Line}
\item \inltt{Spline}
\item \inltt{Spline} (through points)
\item \inltt{BSpline} (i.e. a B\'{e}zier curve)
\item \inltt{Ellipse} (arc): as defined in Gmsh's OpenCASCADE kernel, the first
point defines the start of the arc, the second point the centre and the fourth
point the end. The third point is not used. The start point along with the centre
point form the major axis and the minor axis is then computed so that the end
point falls onto the arc. The major axis must always be greater or equal to the
minor axis.
\item \inltt{Circle} (arc): the circle is a special case of the ellipse where the
third point is skipped. The distances between the start and end points and the
centre point must be equal or an error will be thrown.
\item \inltt{Line Loop}
\item \inltt{Plane Surface}
\end{itemize}
......@@ -950,7 +960,7 @@ At the present time, NekMesh does not support the full scripting capabilities of
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.
splines and arcs.
%%% Local Variables:
......
......@@ -41,6 +41,7 @@ SET(FieldUtilsHeaders
ProcessModules/ProcessPointDataToFld.h
ProcessModules/ProcessPrintFldNorms.h
ProcessModules/ProcessScaleInFld.h
ProcessModules/ProcessStreamFunction.h
ProcessModules/ProcessSurfDistance.h
ProcessModules/ProcessVorticity.h
ProcessModules/ProcessScalGrad.h
......@@ -95,6 +96,7 @@ SET(FieldUtilsSources
ProcessModules/ProcessScaleInFld.cpp
ProcessModules/ProcessVorticity.cpp
ProcessModules/ProcessScalGrad.cpp
ProcessModules/ProcessStreamFunction.cpp
ProcessModules/ProcessSurfDistance.cpp
ProcessModules/ProcessMultiShear.cpp
ProcessModules/ProcessWSS.cpp
......
......@@ -406,6 +406,7 @@ vector<IsoSharedPtr> ProcessIsoContour::ExtractContour(
vector<Array<OneD, int> > ptsConn;
m_f->m_fieldPts->GetConnectivity(ptsConn);
for(int zone = 0; zone < ptsConn.size(); ++zone)
{
IsoSharedPtr iso;
......@@ -529,13 +530,18 @@ vector<IsoSharedPtr> ProcessIsoContour::ExtractContour(
}
}
}
iso->SetNTris(n);
// condense the information in this elemental extraction.
iso->Condense();
returnval.push_back(iso);
if(n)
{
iso->SetNTris(n);
// condense the information in this elemental extraction.
iso->Condense();
returnval.push_back(iso);
}
}
return returnval;
}
......@@ -871,7 +877,7 @@ void Iso::GlobalCondense(vector<IsoSharedPtr> &iso, bool verbose)
for(id1 = 0; id1 < result.size(); ++id1)
{
NekDouble dist = bg::distance(queryPoint, result[id1].first);
if(dist*dist<SQ_PNT_TOL) // same point
if(dist*dist < NekConstants::kNekZeroTol) // same point
{
id2 = result[id1].second;
samept.insert(id2);
......@@ -953,14 +959,14 @@ void Iso::GlobalCondense(vector<IsoSharedPtr> &iso, bool verbose)
bool operator == (const IsoVertex& x, const IsoVertex& y)
{
return ((x.m_x-y.m_x)*(x.m_x-y.m_x) + (x.m_y-y.m_y)*(x.m_y-y.m_y) +
(x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? true:false;
(x.m_z-y.m_z)*(x.m_z-y.m_z) < NekConstants::kNekZeroTol)? true:false;
}
// define != if point is outside 1e-8
bool operator != (const IsoVertex& x, const IsoVertex& y)
{
return ((x.m_x-y.m_x)*(x.m_x-y.m_x) + (x.m_y-y.m_y)*(x.m_y-y.m_y) +
(x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? 0:1;
(x.m_z-y.m_z)*(x.m_z-y.m_z) < NekConstants::kNekZeroTol)? 0:1;
}
void Iso::Smooth(int n_iter, NekDouble lambda, NekDouble mu)
......
......@@ -44,8 +44,6 @@ namespace Nektar
namespace FieldUtils
{
const NekDouble SQ_PNT_TOL=1e-16;
class Iso
{
public:
......@@ -160,13 +158,13 @@ class Iso
m_condensed = false;
m_nvert = 0;
m_fields.resize(nfields);
// set up initial vectors to be 10000 long
m_x.resize(10000);
m_y.resize(10000);
m_z.resize(10000);
// set up initial vectors to be 1000 long
m_x.resize(1000);
m_y.resize(1000);
m_z.resize(1000);
for(int i = 0; i < m_fields.size(); ++i)
{
m_fields[i].resize(10000);
m_fields[i].resize(1000);
}
};
......
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessStreamFunction.cpp
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Computes stream-function for 2D field.
//
////////////////////////////////////////////////////////////////////////////////
#include "ProcessStreamFunction.h"
using namespace std;
namespace Nektar
{
namespace FieldUtils
{
ModuleKey ProcessStreamFunction::className =
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "streamfunction"),
ProcessStreamFunction::create,
"Computes stream-function for 2D field.");
ProcessStreamFunction::ProcessStreamFunction(FieldSharedPtr f)
: ProcessModule(f)
{
m_f->m_declareExpansionAsContField = true;
}
ProcessStreamFunction::~ProcessStreamFunction()
{
}
void ProcessStreamFunction::Process(po::variables_map &vm)
{
int nfields = m_f->m_variables.size();
// Append field name
m_f->m_variables.push_back("StreamFunc");
// Skip in case of empty partition
if (m_f->m_exp[0]->GetNumElmts() == 0)
{
return;
}
// Check if dimension is correct
int expdim = m_f->m_graph->GetMeshDimension();
int spacedim = expdim + m_f->m_numHomogeneousDir;
ASSERTL0(spacedim == 2, "Stream function can only be obtained in 2D.");
// Allocate arrays
int npoints = m_f->m_exp[0]->GetNpoints();
Array<OneD, NekDouble> vx(npoints);
Array<OneD, NekDouble> uy(npoints);
Array<OneD, NekDouble> vortNeg(npoints);
// Resize expansion
m_f->m_exp.resize(nfields + 1);
m_f->m_exp[nfields] = m_f->AppendExpList(m_f->m_numHomogeneousDir);
// Calculate vorticity: -W_z = Uy - Vx
m_f->m_exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
m_f->m_exp[0]->GetPhys(),uy);
m_f->m_exp[1]->PhysDeriv(MultiRegions::DirCartesianMap[0],
m_f->m_exp[1]->GetPhys(),vx);
Vmath::Vsub(npoints, uy, 1, vx, 1, vortNeg, 1);
// Calculate the Stream Function as the solution of the
// Poisson equation: \nabla^2 StreamFunction = -Vorticity
// Use same boundary conditions as v
StdRegions::ConstFactorMap factor;
factor[StdRegions::eFactorLambda] = 0.0;
m_f->m_exp[1]->HelmSolve(vortNeg,
m_f->m_exp[nfields]->UpdateCoeffs(),
NullFlagList, factor);
m_f->m_exp[1]->BwdTrans(m_f->m_exp[nfields]->GetCoeffs(),
m_f->m_exp[nfields]->UpdatePhys());
}
}
}
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessStreamFunction.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Computes stream-function for 2D field.
//
////////////////////////////////////////////////////////////////////////////////
#ifndef FIELDUTILS_PROCESSSTREAMFUNCTION
#define FIELDUTILS_PROCESSSTREAMFUNCTION
#include "../Module.h"
namespace Nektar
{
namespace FieldUtils
{
/**
* @brief This processing module calculates the stream function of a 2D field
* and adds it as an extra-field to the output file.
*/
class ProcessStreamFunction : public ProcessModule
{
public:
/// Creates an instance of this class
static boost::shared_ptr<Module> create(FieldSharedPtr f)
{
return MemoryManager<ProcessStreamFunction>::AllocateSharedPtr(f);
}
static ModuleKey className;
ProcessStreamFunction(FieldSharedPtr f);
virtual ~ProcessStreamFunction();
/// Write mesh to output file.
virtual void Process(po::variables_map &vm);
virtual std::string GetModuleName()
{
return "ProcessStreamFunction";
}
virtual std::string GetModuleDescription()
{
return "Calculating stream-function";
}
virtual ModulePriority GetModulePriority()
{
return eModifyExp;
}
};
}
}
#endif
......@@ -90,6 +90,12 @@ CommMpi::CommMpi(MPI_Comm pComm) : Comm()
*/
CommMpi::~CommMpi()
{
int flag;
MPI_Finalized(&flag);
if (!flag && m_comm != MPI_COMM_WORLD)
{
MPI_Comm_free(&m_comm);
}
}
/**
......
......@@ -2892,7 +2892,6 @@ namespace Nektar
Array<OneD, NekDouble> &bndElmt)
{
int n, cnt, nq;
Array<OneD, NekDouble> tmp1, tmp2;
Array<OneD, int> ElmtID,EdgeID;
GetBoundaryToElmtMap(ElmtID,EdgeID);
......@@ -2920,8 +2919,8 @@ namespace Nektar
{
nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
Vmath::Vcopy(nq, tmp1 = phys + offsetPhys, 1,
tmp2 = bndElmt + offsetElmt, 1);
Vmath::Vcopy(nq, &phys[offsetPhys], 1,
&bndElmt[offsetElmt], 1);
offsetElmt += nq;
}
}
......
......@@ -412,6 +412,9 @@ TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
map<int, string> points;
map<int, string> lines;
map<int, string> splines;
map<int, string> bsplines;
map<int, string> circles;
map<int, string> ellipses;
map<int, string> loops;
map<int, string> surfs;
......@@ -465,6 +468,18 @@ TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
{
splines[id] = var;
}
else if (boost::iequals(type, "BSpline"))
{
bsplines[id] = var;
}
else if (boost::iequals(type, "Circle"))
{
circles[id] = var;
}
else if (boost::iequals(type, "Ellipse"))
{
ellipses[id] = var;
}
else if (boost::iequals(type, "Line Loop"))
{
// line loops sometimes have negative entries for gmsh
......@@ -522,6 +537,109 @@ TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
BRepBuilderAPI_MakeEdge em(curve);
cEdges[it->first] = em.Edge();
}
for (it = bsplines.begin(); it != bsplines.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
TColgp_Array1OfPnt pointArray(0, data.size() - 1);
for (int i = 0; i < data.size(); i++)
{
pointArray.SetValue(i, cPoints[data[i]]);
}
Handle(Geom_BezierCurve) curve = new Geom_BezierCurve(pointArray);
BRepBuilderAPI_MakeEdge em(curve);
cEdges[it->first] = em.Edge();
}
for (it = circles.begin(); it != circles.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
ASSERTL0(data.size() == 3, "Wrong definition of circle arc");
gp_Pnt start = cPoints[data[0]];
gp_Pnt centre = cPoints[data[1]];
gp_Pnt end = cPoints[data[2]];
NekDouble r1 = start.Distance(centre);
NekDouble r2 = end.Distance(centre);
ASSERTL0(fabs(r1 - r2) < 1e-7, "Non-matching radii");
gp_Circ c;
c.SetLocation(centre);
c.SetRadius(r1);
Handle(Geom_Circle) gc = new Geom_Circle(c);
ShapeAnalysis_Curve sac;
NekDouble p1, p2;
sac.Project(gc, start, 1e-8, start, p1);
sac.Project(gc, end, 1e-8, end, p2);
// Make sure the arc is always of length less than pi
if ((p1 > p2) ^ (fabs(p2 - p1) > M_PI))
{
swap(p1, p2);
}
Handle(Geom_TrimmedCurve) tc = new Geom_TrimmedCurve(gc, p1, p2, false);
BRepBuilderAPI_MakeEdge em(tc);
em.Build();
cEdges[it->first] = em.Edge();
}
for (it = ellipses.begin(); it != ellipses.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
ASSERTL0(data.size() == 4, "Wrong definition of ellipse arc");
gp_Pnt start = cPoints[data[0]];
gp_Pnt centre = cPoints[data[1]];
// data[2] useless??
gp_Pnt end = cPoints[data[3]];
NekDouble major = start.Distance(centre);
gp_Vec v1(centre, start);
gp_Vec v2(centre, end);
gp_Vec vx(1.0, 0.0, 0.0);
NekDouble angle = v1.Angle(vx);
// Check for negative rotation
if (v1.Y() < 0)
{
angle *= -1;
}
v2.Rotate(gp_Ax1(), angle);
NekDouble minor = fabs(v2.Y() / sqrt(1.0 - v2.X() * v2.X() / (major * major)));
gp_Elips e;
e.SetLocation(centre);
e.SetMajorRadius(major);
e.SetMinorRadius(minor);