Commit 57f65371 authored by Douglas Serson's avatar Douglas Serson
Browse files

Merge branch 'feature/nodal-cleanup' into 'master'

Rework nodal utilities and add support for nodal prisms

This MR refactors the nodal points utilities that are found in the `NodalUtil.h` and `NodalUtil.cpp` files, and introduces a new `NodalUtil` class which handles the generation of Vandermonde matrices, integration weights, interpolation matrices and derivative matrices for arbitrary sets of points on the triangle, tetrahedron and prism elements. The latter element is newly supported.

This branch:
- Significantly reduces the amount of code from sporadic function calls to a more organised class structure, eliminating most duplicate code
- Fixes a bug in the generation of nodal tetrahedral weights
- Enables the use of nodal prism elements
- Adds a lot of doxygen documentation to describe the new structure
- Replaces the NodalUtil demos with a single NodalDemo, which is now tested in the regression tests

This should be seen as a stepping-stone branch, depending on the outcome of the working being done in the feature/basis-refactor branch, but will hopefully make any further efforts easier.

See merge request !660
parents 27cbe3d9 d847d868
......@@ -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:**
......
......@@ -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()); ++i )
{
for(int j=0; j< int(y.GetRows()); ++j)
{
error(i,j) = LibUtilities::MakeRound(1e15 * fabs(numericXDerivative(i,j) - matVx(i, j)))/1e15;
// Compute relative error
relativeErrorDx(i,j) = (matVx(i,j) - numericXDerivative(i,j))/matVx(i,j);
if( fabs(matVx(i,j)) < numeric_limits<double>::epsilon() )
{
relativeErrorDx(i,j) = matVx(i,j) - numericXDerivative(i,j);
}
}
}
cout << "\n\n\n**************** NodalTetElec X Derivative Floating Point Error Precision ***************" << endl;
cout << "\n Result of NumericXDerivative = \n" << MatrixToString(numericXDerivative) << endl;
cout << "\n Result of exact XDerivative = \n" << MatrixToString(matVx) << endl;
cout << "\n Result of D matrix = \n" << MatrixToString(xderivativeMat) << endl;
cout << "epsilon = \n" << epsilonDx <<endl;
cout << "\n relativeError : X Derivative = \n" << MatrixToString(relativeErrorDx) << endl;
cout << "\n Error : abs(exact - NumericXDerivative) = \n" << MatrixToString(error) << endl;
cout << "\n --------------- End of X Derivative Matrix Demo --------------------------" << endl;
// /////////////////////////////////////////////
// Y Derivative
//
boost::shared_ptr<NekMatrix<NekDouble> > Dyptr = points->GetD(yDir);
const NekMatrix<NekDouble> & yderivativeMat = *Dyptr;
NekMatrix<NekDouble> matVy = GetTetYDerivativeOfMonomialVandermonde(x, y, z);
NekMatrix<NekDouble> numericYDerivative = yderivativeMat * matVnn;
NekMatrix<NekDouble> errorDy(x.GetRows(), y.GetRows());
NekMatrix<NekDouble> relativeErrorDy(x.GetRows(), y.GetRows());
long double epsilonDy = 1e-14;
for(int i=0; i< int(x.GetRows()); ++i )
{
for(int j=0; j< int(y.GetRows()); ++j)
{
error(i,j) = LibUtilities::MakeRound(1e15 * fabs(numericXDerivative(i,j) - matVx(i, j)))/1e15;
// Compute relative error
relativeErrorDy(i,j) = (matVy(i,j) - numericYDerivative(i,j))/matVy(i,j);
if( fabs(matVy(i,j)) < numeric_limits<double>::epsilon() )
{
relativeErrorDy(i,j) = matVy(i,j) - numericYDerivative(i,j);
}
}
}
cout << "\n\n\n*************** NodalTetElec Y Derivative Floating Point Error Precision **************" << endl;
cout << "\n Result of NumericYDerivative = \n" << MatrixToString(numericYDerivative) << endl;
cout << "\n Result of exact Y Derivative = \n" << MatrixToString(matVy) << endl;
cout << "\n Result of D matrix = \n" << MatrixToString(yderivativeMat) << endl;
cout << "epsilon = \n" << epsilonDy <<endl;
cout << "\n relativeError : Y Derivative = \n" << MatrixToString(relativeErrorDy) << endl;
cout << "\n Error : abs(exact - NumericYDerivative) = \n" << MatrixToString(error) << endl;
cout << "\n --------------- End of Y Derivative Matrix --------------------------" << endl;
// /////////////////////////////////////////////
// Z Derivative
//
boost::shared_ptr<NekMatrix<NekDouble> > Dzptr = points->GetD(zDir);
const NekMatrix<NekDouble> & zderivativeMat = *Dzptr;
NekMatrix<NekDouble> matVz = GetTetZDerivativeOfMonomialVandermonde(x, y, z);
NekMatrix<NekDouble> numericZDerivative = zderivativeMat * matVnn;
NekMatrix<NekDouble> errorDz(x.GetRows(), y.GetRows());
NekMatrix<NekDouble> relativeErrorDz(x.GetRows(), y.GetRows());
long double epsilonDz = 1e-14;
for(int i=0; i< int(x.GetRows()); ++i )
{
for(int j=0; j< int(y.GetRows()); ++j)
{
error(i,j) = LibUtilities::MakeRound(1e15 * fabs(numericZDerivative(i,j) - matVz(i, j)))/1e15;
// Compute relative error