Commit 3d1ff98a authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'master' into feature/qmetric_scaled

parents 3dbf72ee f238a4b9
......@@ -23,6 +23,7 @@ v4.4.0
to boundary conditions (!615)
- Allow expansions to be loaded directly from field file (!617)
- New options for load balancing (DOF or BOUNDARY) in mesh partitioner (!617)
- Rework nodal utilities to support nodal prismatic elements (!660)
- Update Body/Field forces at each timestep (!665)
**ADRSolver:**
......@@ -46,8 +47,14 @@ v4.4.0
- Add STL surface writer module (!668)
- New module for inserting an alternate high-order surface into the working
mesh (!669)
- Add curve projection routines to CAD system (!697)
- Improvements to mesh linearisation module (!659)
- Add support for Gmsh high-order output (!679)
- Move CAD classes to factory format (!676)
- Add module to check topology of the mesh along with boundary connectivity
to detect problems such as hanging nodes (!691)
- Add option to `linearise` module to linearise only prisms (!688)
- Add option to `linearise` to use element quality (!690)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
......@@ -15,7 +15,7 @@ IF(NEKTAR_USE_MESHGEN)
ELSE()
SET(BUILD_OCC ON)
ENDIF()
OPTION(THIRDPARTY_BUILD_OCC "Build OpenCascade library from ThirdParty."
${BUILD_OCC})
......@@ -33,7 +33,7 @@ IF(NEKTAR_USE_MESHGEN)
IF(WIN32)
MESSAGE(SEND_ERROR "Cannot currently use OpenCascade with Nektar++ on Windows")
ENDIF()
EXTERNALPROJECT_ADD(
opencascade-6.9
PREFIX ${TPSRC}
......
......@@ -542,17 +542,32 @@ The module parameters are:
\subsection{Linearisation}
The ability to remove all the high-order information in a mesh can be useful
at times.
The ability to remove all the high-order information in a mesh can be useful at
times when mesh generation is tricky or impossible in the presence of curvature
To do this in NekMesh use the command:
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m linearise high-order-mesh.xml linear-mesh.xml
NekMesh -m linearise:all high-order-mesh.xml linear-mesh.xml
\end{lstlisting}
The output will contain only the linear mesh information, all curved information
is removed.
is removed. Alternatively
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m linearise:invalid high-order-mesh.xml linear-mesh.xml
\end{lstlisting}
attempts to remove curvature from elements only where necessary. This is a
simple algorithm that removes curvature from invalid elements and repeats until
all elements are valid. Either \inlsh{all} or \inlsh{invalid} must be specified.
\begin{itemize}
\item \inlsh{all}: remove curvature from all elements.
\item \inlsh{invalid}: remove curvature from invalid elements.
\item \inlsh{prismonly}: consider only prisms when removing curvature. This is
useful in the presence of a prismatic boundary layer.
\end{itemize}
\subsection{Extracting interface between tetrahedra and prismatic elements}
......@@ -603,6 +618,19 @@ may issue the command
module in combination with the splitting module described earlier.
\end{notebox}
\subsection{Link Checking}
It is quite possible that a mesh contains some sort of hanging entity or
element connectivity error. The check link module is a fast check that, a)
elements are correctly connected and b) the boundary entities (composites)
match the interior domain:
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m linkcheck mesh.xml mesh2.xml
\end{lstlisting}
This module should be added to the module chain for any complex if the user
suspects there may be a mesh issue. The module will print a warning if there is.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......
......@@ -11,20 +11,7 @@ SET(PartitionAnalyseSources
SET(FoundationSources
FoundationDemo.cpp)
SET(NodalTriFeketeSources
NodalTriFeketeDemo.cpp)
SET(NodalTriElecSources
NodalTriElecDemo.cpp)
SET(NodalTriEvenlySpacedSources
NodalTriEvenlySpacedDemo.cpp)
SET(NodalTetEvenlySpacedSources
NodalTetEvenlySpacedDemo.cpp)
SET(NodalTetElecSources
NodalTetElecDemo.cpp)
SET(NodalDemoSources NodalDemo.cpp)
SET(TimeIntegrationDemoSources
TimeIntegrationDemo.cpp)
......@@ -43,24 +30,22 @@ TARGET_LINK_LIBRARIES(PartitionAnalyse LibUtilities)
ADD_NEKTAR_EXECUTABLE(FoundationDemo demos FoundationSources )
TARGET_LINK_LIBRARIES(FoundationDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTriFeketeDemo demos NodalTriFeketeSources )
TARGET_LINK_LIBRARIES(NodalTriFeketeDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTriElecDemo demos NodalTriElecSources)
TARGET_LINK_LIBRARIES(NodalTriElecDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTriEvenlySpacedDemo demos NodalTriEvenlySpacedSources)
TARGET_LINK_LIBRARIES(NodalTriEvenlySpacedDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTetEvenlySpacedDemo demos NodalTetEvenlySpacedSources)
TARGET_LINK_LIBRARIES(NodalTetEvenlySpacedDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalTetElecDemo demos NodalTetElecSources)
TARGET_LINK_LIBRARIES(NodalTetElecDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(NodalDemo demos NodalDemoSources)
TARGET_LINK_LIBRARIES(NodalDemo LibUtilities)
ADD_NEKTAR_EXECUTABLE(TimeIntegrationDemo demos TimeIntegrationDemoSources)
TARGET_LINK_LIBRARIES(TimeIntegrationDemo LibUtilities)
ADD_NEKTAR_TEST(NodalDemo_Tri_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tri_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tri_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Prism_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Prism_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Prism_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Tet_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tet_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tet_Interp_P7)
IF(NEKTAR_USE_HDF5)
ADD_NEKTAR_EXECUTABLE(FieldIOBenchmarker demos FieldIOBenchmarkerSources)
TARGET_LINK_LIBRARIES(FieldIOBenchmarker LibUtilities)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: NodalDemo.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: Demo for testing functionality of NodalUtil classes
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <string>
#include <boost/program_options.hpp>
#include <boost/algorithm/string.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/NodalUtil.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/BasicUtils/ShapeType.hpp>
using namespace std;
using namespace Nektar;
using namespace Nektar::LibUtilities;
namespace po = boost::program_options;
int main(int argc, char *argv[])
{
po::options_description desc("Available options");
desc.add_options()
("help,h", "Produce this help message.")
("integral,i", "Evaluate the integral of a known function and "
"return error.")
("order,o", po::value<int>(),
"Order to evaluate at.")
("type,t", po::value<int>(),
"Points type of the nodal element.")
("deriv,d", "Evaluate the derivative of a known function and "
"return error.")
("interp,p", po::value<string>(),
"Interpolate function at a known point of the form"
"x,y,z and return error.");
po::variables_map vm;
try
{
po::store(po::command_line_parser(argc, argv).
options(desc).run(), vm);
po::notify(vm);
}
catch (const std::exception& e)
{
cerr << e.what() << endl;
cerr << desc;
return 1;
}
if (!vm.count("order"))
{
cerr << "Must supply an element order." << endl;
return 1;
}
int order = vm["order"].as<int>();
map<PointsType, ShapeType> nodalTypes;
map<PointsType, ShapeType>::iterator nIt;
nodalTypes[eNodalTriElec] = eTriangle;
nodalTypes[eNodalTriFekete] = eTriangle;
nodalTypes[eNodalTriEvenlySpaced] = eTriangle;
nodalTypes[eNodalTetEvenlySpaced] = eTetrahedron;
nodalTypes[eNodalTetElec] = eTetrahedron;
nodalTypes[eNodalPrismEvenlySpaced] = ePrism;
if (!vm.count("type"))
{
cerr << "Must supply a points type. Valid points types are:" << endl;
for (nIt = nodalTypes.begin(); nIt != nodalTypes.end(); ++nIt)
{
cerr << " - " << (int)nIt->first << " ("
<< kPointsTypeStr[nIt->first] << ")" << endl;
}
return 1;
}
PointsType pointsType = (PointsType)vm["type"].as<int>();
if (nodalTypes.find(pointsType) == nodalTypes.end())
{
cerr << "Invalid points type. Valid points types are:" << endl;
for (nIt = nodalTypes.begin(); nIt != nodalTypes.end(); ++nIt)
{
cerr << " - " << (int)nIt->first << " ("
<< kPointsTypeStr[nIt->first] << ")" << endl;
}
return 1;
}
// Generate nodal points.
PointsKey pointsKey(order + 1, pointsType);
PointsSharedPtr points = PointsManager()[pointsKey];
Array<OneD, NekDouble> r, s, t;
points->GetPoints(r, s, t);
// Generate nodal utility.
ShapeType shape = nodalTypes[pointsType];
NodalUtil *util = NULL;
if (shape == eTriangle)
{
util = new NodalUtilTriangle(order, r, s);
}
else if (shape == eTetrahedron)
{
util = new NodalUtilTetrahedron(order, r, s, t);
}
else if (shape == ePrism)
{
util = new NodalUtilPrism(order, r, s, t);
}
ASSERTL1(util, "Unknown shape type!");
const int nPoints = r.num_elements();
const int dim = shape == eTriangle ? 2 : 3;
if (vm.count("integral"))
{
NekVector<NekDouble> weights = util->GetWeights();
NekDouble integral = 0.0;
for (int i = 0; i < nPoints; ++i)
{
NekDouble integrand = dim == 2 ?
exp(r[i]) * sin(s[i]) : exp(r[i] + s[i] + t[i]);
integral += weights[i] * integrand;
}
NekDouble exact = 0.0;
switch(shape)
{
case eTriangle:
exact = -0.5 * (sin(1.0) + cos(1.0) + M_E * M_E *
(sin(1.0) - cos(1.0))) / M_E;
break;
case eTetrahedron:
exact = 1.0 / M_E - 1.0 / M_E / M_E / M_E;
break;
case ePrism:
exact = M_E - 1.0 / M_E / M_E / M_E;
break;
default:
exact = 0.0;
break;
}
cout << "L inf error : " << fabs(exact - integral) << endl;
}
else if (vm.count("deriv"))
{
Array<OneD, NekVector<NekDouble> > exact(dim);
NekVector<NekDouble> input(nPoints), output(nPoints);
for (int i = 0; i < dim; ++i)
{
exact[i] = NekVector<NekDouble>(nPoints);
}
if (dim == 2)
{
// Exact solution: u(x,y) = sin(x) * cos(y)
for (int i = 0; i < nPoints; ++i)
{
input[i] = sin(r[i]) * cos(s[i]);
exact[0][i] = cos(r[i]) * cos(s[i]);
exact[1][i] = -sin(r[i]) * sin(s[i]);
}
}
else
{
// Exact solution: u(x,y) = sin(x) * cos(y) * sin(z)
for (int i = 0; i < nPoints; ++i)
{
input[i] = sin(r[i]) * cos(s[i]) * sin(t[i]);
exact[0][i] = cos(r[i]) * cos(s[i]) * sin(t[i]);
exact[1][i] = -sin(r[i]) * sin(s[i]) * sin(t[i]);
exact[2][i] = sin(r[i]) * cos(s[i]) * cos(t[i]);
}
}
char *vars[3] = { "u", "v", "w" };
for (int i = 0; i < dim; ++i)
{
boost::shared_ptr<NekMatrix<NekDouble> > deriv =
util->GetDerivMatrix(i);
output = *deriv * input;
NekVector<NekDouble> tmp = output - exact[i];
cout << "L 2 error (variable " << vars[i] << ") : " << tmp.L2Norm()
<< endl;
}
}
else if (vm.count("interp"))
{
string pointStr = vm["interp"].as<string>();
vector<string> point;
boost::split(point, pointStr, boost::is_any_of(","));
if (point.size() != dim)
{
cerr << "Point " << pointStr << " does not have correct dimensions"
<< endl;
return 1;
}
Array<OneD, Array<OneD, NekDouble> > tmp(dim);
for (int i = 0; i < dim; ++i)
{
tmp[i] = Array<OneD, NekDouble>(1);
try
{
tmp[i][0] = boost::lexical_cast<NekDouble>(point[i]);
}
catch(boost::bad_lexical_cast &)
{
cerr << "Could not convert " << point[i] << " to a coordinate"
<< endl;
return 1;
}
}
NekVector<NekDouble> input(nPoints);
for (int i = 0; i < nPoints; ++i)
{
input[i] = dim == 2 ? exp(r[i]) * exp(s[i]) :
exp(r[i]) * exp(s[i]) * exp(t[i]);
}
boost::shared_ptr<NekMatrix<NekDouble> > interp =
util->GetInterpolationMatrix(tmp);
NekVector<NekDouble> output = *interp * input;
NekDouble exact = dim == 2 ? exp(tmp[0][0]) * exp(tmp[1][0]) :
exp(tmp[0][0]) * exp(tmp[1][0]) * exp(tmp[2][0]);
cout << "L inf error : " << fabs(exact - output[0]) << endl;
}
else
{
cerr << "You must specify one of --integral, --deriv or --interp"
<< endl;
cerr << desc;
return 1;
}
return 0;
}
#include <iostream>
#include <iomanip>
using namespace std;
#include <LibUtilities/Foundations/NodalUtil.h>
#include <LibUtilities/Foundations/NodalTetElec.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Points.h>
using namespace Nektar;
using namespace boost;
using namespace Nektar::LibUtilities;
int main(int argc, char *argv[]){
// Argument check: Display a help message if the count is wrong
if(argc != 3)
{
cerr << "Usage: NodalTetElecDemo Points3D-Type nPtsPerSide" << endl;
cerr << "Where type is an integer value which dictates the basis as:" << endl;
for(int i=0; i<SIZE_PointsType; ++i)
{
cerr << setw(30) << kPointsTypeStr[i] << " =" << i << endl;
}
cerr << "Note type = 1 ~ 14 are for one dimensional basis, and 15 and 16 are for two dimensional basis " << endl;
cerr << "\n Example: NodalTetElecDemo 17 4" << endl;
cerr << "\n\t Tests 3D NodalTetElec on 4 points per side" << endl;
return 1; // Aborts main() function
}
// Read in the type for the points from the caller
PointsType pointsType = (PointsType) atoi(argv[1]);
if(pointsType == eNoPointsType)
{
cerr << "pointsType = " << pointsType << endl;
cerr << "kPointsTypeStr[" <<pointsType<< "]=" << kPointsTypeStr[pointsType] << endl;
ErrorUtil::Error(ErrorUtil::efatal,__FILE__, __LINE__,
"No Points Type requested",0);
}
// Read in the number of points per side at which the function will be evaluated
int nPtsPerSide = atoi(argv[2]);
// Show the set up to the user
cout << "Points Type: " << kPointsTypeStr[pointsType] << endl;
cout << "Number of points per side: " << nPtsPerSide << endl;
// Display the example test function to the user
if( pointsType == eNodalTetElec)
{
cout << "Uniform grid on a Tetrahedron points" << endl;
}
int degree = nPtsPerSide-1;
boost::shared_ptr<Points<NekDouble> > points = PointsManager()[PointsKey(nPtsPerSide, pointsType)];
NodalTetElec *nTete = dynamic_cast<NodalTetElec*>(points.get());
Array<OneD, const NekDouble> ax, ay, az;
nTete->GetPoints(ax,ay, az);
NekVector<NekDouble> x = ToVector(ax), y = ToVector(ay), z = ToVector(az);
// /////////////////////////////////////////////
// Test Interpolation
//
// Make a uniform grid on the Tetrahedron
int nGridPtsPerSide =6;
int nGridPts = nGridPtsPerSide * (nGridPtsPerSide + 1) * (nGridPtsPerSide + 2) / 6;
Array<OneD, NekDouble> axi(nGridPts), ayi(nGridPts), azi(nGridPts);
for( int i = 0, n = 0; i < nGridPtsPerSide; ++i )
{
for( int j = 0; j < nGridPtsPerSide - i; ++j )
{
for(int k = 0; k < nGridPtsPerSide - i - j; ++k, ++n )
{
axi[n] = -1.0 + 2.0*j / (nGridPtsPerSide - 1);
ayi[n] = -1.0 + 2.0*i / (nGridPtsPerSide - 1);
azi[n] = -1.0 + 2.0*k / (nGridPtsPerSide - 1);
}
}
}
// Done with constructing the uniform grid of points on the Tetrahedron
NekVector<NekDouble> xi = ToVector(axi), yi = ToVector(ayi), zi = ToVector(azi);
boost::shared_ptr<NekMatrix<NekDouble> > Iptr = nTete->GetI(axi, ayi, azi);
const NekMatrix<NekDouble> & interpMat = *Iptr;
NekMatrix<NekDouble> matVnn = GetMonomialVandermonde(x, y, z);
NekMatrix<NekDouble> matVmn = GetMonomialVandermonde(xi, yi, zi, degree);
NekMatrix<NekDouble> matVmnTilde = interpMat * matVnn;
NekMatrix<NekDouble> err(xi.GetRows(), x.GetRows());
NekMatrix<NekDouble> relativeError(GetSize(xi), GetSize(x));
long double epsilon = 1e-15;
for( int i = 0; i < int(xi.GetRows()); ++i )
{
for( int j = 0; j < int(x.GetRows()); ++j )
{
err(i,j) = LibUtilities::MakeRound(1e16 * fabs(matVmnTilde(i,j) - matVmn(i,j)))/1e16;
// Compute relative error
relativeError(i,j) = (matVmn(i,j) - matVmnTilde(i,j))/matVmn(i,j);
if( fabs(matVmn(i,j)) < numeric_limits<double>::epsilon() )
{
relativeError(i,j) = matVmn(i,j) - matVmnTilde(i,j);
}
}
}
cout << "\n\n\n************************* NodalTetElec Interpolation Demo *******************************" << endl;
cout << "\n Result of NumericInterpolation = \n" << MatrixToString(matVmnTilde) << endl;
cout << "\n Result of Exact Interpolation = \n" << MatrixToString(matVmn) << endl;
cout << "\n Result of I matrix = \n" << MatrixToString(interpMat) << endl;
cout << "\n epsilon = \n" << epsilon << endl;
cout << "\n relativeError : Interpolation = \n" << MatrixToString(relativeError) << endl;
cout << "\n Error : abs(NumericInterpolation - exact) = \n" << MatrixToString(err) << endl;
cout << "---------------------------------- End of Interpolation Demo ------------------------------------" << endl;
// /////////////////////////////////////////////
// X Derivative
//
boost::shared_ptr<NekMatrix<NekDouble> > Dxptr = points->GetD(xDir);
const NekMatrix<NekDouble> & xderivativeMat = *Dxptr;
NekMatrix<NekDouble> matVx = GetTetXDerivativeOfMonomialVandermonde(x, y, z);
NekMatrix<NekDouble> numericXDerivative = xderivativeMat * matVnn;
NekMatrix<NekDouble> error(x.GetRows(), y.GetRows());
NekMatrix<NekDouble> relativeErrorDx(x.GetRows(), y.GetRows());
long double epsilonDx = 1e-14;
for(int i=0; i< int(x.GetRows