Skip to content

GitLab

Nektar
Nektar
Commit edd2ac0a authored by Chris Cantwell's avatar Chris Cantwell
Browse files
Options

Merge branch 'master' into feature/scotch 

Conflicts:
	utilities/PreProcessing/MeshConvert/CMakeLists.txt
parents 2e95e5d3 742629f3
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 82256 additions and 358 deletions
docs/html/doxygen/doxygen.bib 0 → 100644
View file @ edd2ac0a
@article{CaYaKiPeSh13,
title={High-order spectral/< i> hp</i> element discretisation for reaction-diffusion problems on surfaces: Application to cardiac electrophysiology},
author={Cantwell, CD and Yakovlev, S and Kirby, RM and Peters, NS and Sherwin, SJ},
journal={Journal of Computational Physics},
year={2013},
publisher={Elsevier}
}
library/LibUtilities/BasicUtils/FieldIO.cpp
View file @ edd2ac0a
......@@ -58,6 +58,7 @@
namespace ptime = boost::posix_time;
namespace ip = boost::asio::ip;
namespace berrc = boost::system::errc;
namespace Nektar
{
......@@ -149,8 +150,6 @@ namespace Nektar
TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
doc.LinkEndChild(decl);
cout << "Writing outfile: " << filename << endl;
TiXmlElement * root = new TiXmlElement("NEKTAR");
doc.LinkEndChild(root);
......@@ -477,8 +476,6 @@ namespace Nektar
ASSERTL0(fileNames.size() == elementList.size(),"Outfile names and list of elements ids does not match");
cout << "Writing multi-file data: " << outFile << endl;
TiXmlElement * root = new TiXmlElement("NEKTAR");
doc.LinkEndChild(root);
......@@ -1084,40 +1081,14 @@ namespace Nektar
fs::path specPath (outname);
// Remove any existing file which is in the way
int existCheck = fs::exists(specPath) ? 1 : 0;
m_comm->AllReduce(existCheck, ReduceMax);
if (existCheck)
try
{
// First remove all files on the root process.
if (m_comm->GetRank() == 0)
{
fs::remove_all(specPath);
}
m_comm->Block();
// Check to see if the files still exist on non-root processes.
int existCheck = rank > 0 && fs::exists(specPath) ? 1 : 0;
m_comm->AllReduce(existCheck, ReduceMax);
// If they do (i.e. a non-shared filesystem) then go through all
// other processors and try to perform removal there. Note this
// could be made quicker by the use of a per-node MPI
// communicator.
if (existCheck > 0)
{
for (int i = 1; i < nprocs; ++i)
{
m_comm->Block();
if (rank == i && fs::exists(specPath))
{
// Recursively remove directories
fs::remove_all(specPath);
}
}
}
fs::remove_all(specPath);
}
catch (fs::filesystem_error& e)
{
ASSERTL0(e.code().value() == berrc::no_such_file_or_directory,
"Filesystem error: " + string(e.what()));
}
// serial processing just add ending.
......@@ -1140,6 +1111,16 @@ namespace Nektar
}
m_comm->AllReduce(elmtnums,LibUtilities::ReduceMax);
// Create the destination directory
try
{
fs::create_directory(specPath);
}
catch (fs::filesystem_error& e)
{
ASSERTL0(false, "Filesystem error: " + string(e.what()));
}
// Collate per-process element lists on root process to generate
// the info file.
if (rank == 0)
......@@ -1164,11 +1145,11 @@ namespace Nektar
filenames.push_back(pad.str());
}
// Create the destination directory
fs::create_directory(specPath);
// Write the Info.xml file
string infofile = LibUtilities::PortablePath(specPath / fs::path("Info.xml"));
string infofile = LibUtilities::PortablePath(
specPath / fs::path("Info.xml"));
cout << "Writing: " << specPath << endl;
WriteMultiFldFileIDs(infofile, filenames, ElementIDs,
fieldmetadatamap);
}
......@@ -1178,18 +1159,6 @@ namespace Nektar
m_comm->Send(0, idlist);
}
// Ensure all processors are aligned for the write and ensure
// target directory has been created by the root process
m_comm->Block();
// Sit in a loop and make sure target directory has been created
int created = 0;
do
{
created = fs::is_directory(specPath) ? 1 : 0;
m_comm->AllReduce(created, ReduceMin);
} while (!created);
// Pad rank to 8char filenames, e.g. P0000000.fld
boost::format pad("P%1$07d.fld");
pad % m_comm->GetRank();
......
library/LibUtilities/BasicUtils/MeshPartition.cpp
View file @ edd2ac0a
......@@ -105,7 +105,6 @@ namespace Nektar
TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
vNew.LinkEndChild(decl);
TiXmlElement* vElmtNektar;
vElmtNektar = new TiXmlElement("NEKTAR");
......@@ -114,7 +113,6 @@ namespace Nektar
vNew.LinkEndChild(vElmtNektar);
std::string dirname = pSession->GetSessionName() + "_xml";
fs::path pdirname(dirname);
......@@ -122,14 +120,10 @@ namespace Nektar
pad % rank;
fs::path pFilename(pad.str());
if(rank == 0)
if(!fs::is_directory(dirname))
{
if(!fs::is_directory(dirname))
{
fs::create_directory(dirname);
}
fs::create_directory(dirname);
}
m_comm->Block();
fs::path fullpath = pdirname / pFilename;
vNew.SaveFile(PortablePath(fullpath));
......@@ -153,8 +147,9 @@ namespace Nektar
std::string dirname = pSession->GetSessionName() + "_xml";
fs::path pdirname(dirname);
std::string vFilename = "P" + boost::lexical_cast<std::string>(i) + ".xml";
fs::path pFilename(vFilename);
boost::format pad("P%1$07d.xml");
pad % i;
fs::path pFilename(pad.str());
fs::path fullpath = pdirname / pFilename;
......
library/LibUtilities/BasicUtils/Tau.hpp 0 → 100644
View file @ edd2ac0a
///////////////////////////////////////////////////////////////////////////////
//
// File Tau.hpp
//
// 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: wrapper of functions around TAU routines
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_LIB_UTILITIES_BASICUTILS_TAU_HPP
#define NEKTAR_LIB_UTILITIES_BASICUTILS_TAU_HPP
namespace Tau
{
extern "C"
{
void Tau_start(const char *name);
void Tau_stop(const char *name);
}
Start(std::string name)
{
Tau_start(name.c_str());
}
Stop(std::string name)
{
Tau_stop(name.c_str());
}
}
#endif //NEKTAR_LIB_UTILITIES_BASICUTILS_TAU_HPP
library/LibUtilities/Foundations/BLPoints.cpp
View file @ edd2ac0a
......@@ -54,6 +54,14 @@ namespace Nektar
"Must set factor in BLPoints key");
if (fabs(r-1.0) < 1e-6)
{
NekDouble tmp = 2.0/(npts-1.0);
for (unsigned int i = 0; i < npts; ++i)
{
m_points[0][i] = -1.0 + i * tmp;
}
}
else
{
NekDouble a = 2.0 * (1.0-r) / (1.0 - pow(r,(double)npts));
m_points[0][0] = -1.0;
......@@ -65,14 +73,6 @@ namespace Nektar
m_points[0][npts-1] = 1.0;
}
else
{
NekDouble tmp = 2.0/(npts-1.0);
for (unsigned int i = 0; i < npts; ++i)
{
m_points[0][i] = -1.0 + i * tmp;
}
}
if (m_pointsKey.GetPointsType() == eBoundaryLayerPointsRev)
{
......
library/LocalRegions/HexExp.cpp
View file @ edd2ac0a
......@@ -1598,8 +1598,10 @@ namespace Nektar
break;
case StdRegions::eLaplacian:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
mkey.GetNVarCoeff())
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
mkey.GetNVarCoeff()||
mkey.ConstFactorExists(
StdRegions::eFactorSVVCutoffRatio))
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
......
library/LocalRegions/PrismExp.cpp
View file @ edd2ac0a
......@@ -1341,7 +1341,10 @@ namespace Nektar
case StdRegions::eLaplacian:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
mkey.GetNVarCoeff() > 0 ||
mkey.ConstFactorExists(
StdRegions::eFactorSVVCutoffRatio))
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
......
library/LocalRegions/PyrExp.cpp
View file @ edd2ac0a
......@@ -451,7 +451,10 @@ namespace Nektar
break;
case StdRegions::eLaplacian:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
mkey.GetNVarCoeff() > 0 ||
mkey.ConstFactorExists(
StdRegions::eFactorSVVCutoffRatio))
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
......
library/MultiRegions/AssemblyMap/AssemblyMapCG3D.cpp
View file @ edd2ac0a
......@@ -1502,7 +1502,7 @@ namespace Nektar
if (iter != extraDirVerts.end() &&
foundExtraVerts.count(meshVertId) == 0)
{
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
bndCondFaceExp->GetVertexMap(k);
int gid = graphVertOffset[
vertReorderedGraphVertId[meshVertId]];
......@@ -1574,11 +1574,11 @@ namespace Nektar
{
for(l = 0; l < nEdgeInteriorCoeffs; ++l)
{
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
edgeInteriorMap[l];
int gid = graphVertOffset[
edgeReorderedGraphVertId[meshEdgeId]]+l;
ExtraDirDof t(loc, gid, 1.0);
ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
m_extraDirDofs[i].push_back(t);
}
foundExtraEdges.insert(meshEdgeId);
......@@ -1655,7 +1655,7 @@ namespace Nektar
{
for (i = 0; i < Tit->second.size(); ++i)
{
boost::get<2>(Tit->second.at(i)) = 1.0/valence[Tit->second.at(i).get<1>()];
boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
}
}
......
library/MultiRegions/DisContField3D.cpp
View file @ edd2ac0a
......@@ -2418,13 +2418,45 @@
if(m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eDirichlet)
{
LibUtilities::Equation condition = boost::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition >(m_bndConditions[i])->m_dirichletCondition;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
string filebcs = boost::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition>(
m_bndConditions[i])->m_filename;
if(filebcs != "")
{
string var = filebcs.substr(
0, filebcs.find_last_of("."));
int len = var.length();
var = var.substr(len-1,len);
cout << "Boundary condition from file:"
<< filebcs << endl;
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
Import(filebcs,FieldDef, FieldData);
// copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[0], locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else
{
LibUtilities::Equation condition = boost::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition >(m_bndConditions[i])->m_dirichletCondition;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
}
else if(m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eNeumann)
......
library/SolverUtils/Filters/FilterAeroForces.cpp
View file @ edd2ac0a
......@@ -175,9 +175,27 @@ namespace Nektar
{
// Open output stream
m_outputStream.open(m_outputFile.c_str());
m_outputStream << "# Time\t Fx (press)\t Fx (visc)\t Fx (tot)\t"
" Fy (press)\t Fy (visc)\t Fy (tot)\t"
" Fz (press)\t Fz (visc)\t Fz (tot)\t";
m_outputStream << "#";
m_outputStream.width(7);
m_outputStream << "Time";
m_outputStream.width(25);
m_outputStream << "Fx (press)";
m_outputStream.width(25);
m_outputStream << "Fx (visc)";
m_outputStream.width(25);
m_outputStream << "Fx (tot)";
m_outputStream.width(25);
m_outputStream << "Fy (press)";
m_outputStream.width(25);
m_outputStream << "Fy (visc)";
m_outputStream.width(25);
m_outputStream << "Fy (tot)";
m_outputStream.width(25);
m_outputStream << "Fz (press)";
m_outputStream.width(25);
m_outputStream << "Fz (visc)";
m_outputStream.width(25);
m_outputStream << "Fz (tot)";
m_outputStream << endl;
}
......@@ -280,9 +298,9 @@ namespace Nektar
// points for each element (hence local).
for(int j = 0; j < dim; ++j)
{
gradU[j] = Array<OneD, NekDouble>(nq);
gradV[j] = Array<OneD, NekDouble>(nq);
gradW[j] = Array<OneD, NekDouble>(nq);
gradU[j] = Array<OneD, NekDouble>(nq,0.0);
gradV[j] = Array<OneD, NekDouble>(nq,0.0);
gradW[j] = Array<OneD, NekDouble>(nq,0.0);
}
// identify boundary of element
......@@ -305,20 +323,20 @@ namespace Nektar
int nbc = bc->GetTotPoints();
// several vectors for computing the forces
Array<OneD, NekDouble> Pb(nbc);
Array<OneD, NekDouble> Pb(nbc,0.0);
for(int j = 0; j < dim; ++j)
{
fgradU[j] = Array<OneD, NekDouble>(nbc);
fgradV[j] = Array<OneD, NekDouble>(nbc);
fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
}
Array<OneD, NekDouble> drag_t(nbc);
Array<OneD, NekDouble> lift_t(nbc);
Array<OneD, NekDouble> drag_p(nbc);
Array<OneD, NekDouble> lift_p(nbc);
Array<OneD, NekDouble> temp(nbc);
Array<OneD, NekDouble> temp2(nbc);
Array<OneD, NekDouble> drag_t(nbc,0.0);
Array<OneD, NekDouble> lift_t(nbc,0.0);
Array<OneD, NekDouble> drag_p(nbc,0.0);
Array<OneD, NekDouble> lift_p(nbc,0.0);
Array<OneD, NekDouble> temp(nbc,0.0);
Array<OneD, NekDouble> temp2(nbc,0.0);
// identify boundary of element .
boundary = BoundarytoTraceID[cnt];
......@@ -442,9 +460,9 @@ namespace Nektar
// points for each element (hence local).
for(int j = 0; j < dim; ++j)
{
gradU[j] = Array<OneD, NekDouble>(nq);
gradV[j] = Array<OneD, NekDouble>(nq);
gradW[j] = Array<OneD, NekDouble>(nq);
gradU[j] = Array<OneD, NekDouble>(nq,0.0);
gradV[j] = Array<OneD, NekDouble>(nq,0.0);
gradW[j] = Array<OneD, NekDouble>(nq,0.0);
}
//identify boundary of element
......@@ -469,24 +487,24 @@ namespace Nektar
int nbc = bc->GetTotPoints();
//several vectors for computing the forces
Array<OneD, NekDouble> Pb(nbc);
Array<OneD, NekDouble> Pb(nbc,0.0);
for(int j = 0; j < dim; ++j)
{
fgradU[j] = Array<OneD, NekDouble>(nbc);
fgradV[j] = Array<OneD, NekDouble>(nbc);
fgradW[j] = Array<OneD, NekDouble>(nbc);
fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
fgradW[j] = Array<OneD, NekDouble>(nbc,0.0);
}
Array<OneD, NekDouble> drag_t(nbc);
Array<OneD, NekDouble> lift_t(nbc);
Array<OneD, NekDouble> side_t(nbc);
Array<OneD, NekDouble> drag_p(nbc);
Array<OneD, NekDouble> lift_p(nbc);
Array<OneD, NekDouble> side_p(nbc);
Array<OneD, NekDouble> temp(nbc);
Array<OneD, NekDouble> temp2(nbc);
Array<OneD, NekDouble> drag_t(nbc,0.0);
Array<OneD, NekDouble> lift_t(nbc,0.0);
Array<OneD, NekDouble> side_t(nbc,0.0);
Array<OneD, NekDouble> drag_p(nbc,0.0);
Array<OneD, NekDouble> lift_p(nbc,0.0);
Array<OneD, NekDouble> side_p(nbc,0.0);
Array<OneD, NekDouble> temp(nbc,0.0);
Array<OneD, NekDouble> temp2(nbc,0.0);
// identify boundary of element .
boundary = BoundarytoTraceID[cnt];
......@@ -638,8 +656,8 @@ namespace Nektar
for(int j = 0; j < dim; ++j)
{
gradU[j] = Array<OneD, NekDouble>(nq);
gradV[j] = Array<OneD, NekDouble>(nq);
gradU[j] = Array<OneD, NekDouble>(nq,0.0);
gradV[j] = Array<OneD, NekDouble>(nq,0.0);
}
boundary = BoundarytoTraceID[cnt];
......@@ -655,13 +673,13 @@ namespace Nektar
::Expansion1D> (BndExp[n]->GetExp(i));
int nbc = bc->GetTotPoints();
Array<OneD, NekDouble> Pb(nbc);
Array<OneD, NekDouble> Pb(nbc,0.0);
Array<OneD, NekDouble> drag_t(nbc);
Array<OneD, NekDouble> lift_t(nbc);
Array<OneD, NekDouble> drag_p(nbc);
Array<OneD, NekDouble> lift_p(nbc);
Array<OneD, NekDouble> temp(nbc);
Array<OneD, NekDouble> drag_t(nbc,0.0);
Array<OneD, NekDouble> lift_t(nbc,0.0);
Array<OneD, NekDouble> drag_p(nbc,0.0);
Array<OneD, NekDouble> lift_p(nbc,0.0);
Array<OneD, NekDouble> temp(nbc,0.0);
boundary = BoundarytoTraceID[cnt];
......@@ -669,8 +687,8 @@ namespace Nektar
for(int j = 0; j < dim; ++j)
{
fgradU[j] = Array<OneD, NekDouble>(nbc);
fgradV[j] = Array<OneD, NekDouble>(nbc);
fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
}
......@@ -746,7 +764,7 @@ namespace Nektar
m_outputStream.width(25);
m_outputStream << setprecision(8) << Fxv;
m_outputStream.width(25);
m_outputStream << setprecision(16) << Fx;
m_outputStream << setprecision(8) << Fx;
m_outputStream.width(25);
m_outputStream << setprecision(8) << Fyp;
......
library/StdRegions/StdExpansion3D.cpp
View file @ edd2ac0a
......@@ -197,7 +197,8 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
if ( mkey.GetNVarCoeff() == 0 &&
!mkey.ConstFactorExists(eFactorSVVCutoffRatio))
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
......
library/StdRegions/StdHexExp.cpp
View file @ edd2ac0a
......@@ -2394,132 +2394,7 @@ namespace Nektar
StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
}
void StdHexExp::v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdMatrixKey &mkey)
{
ASSERTL0(false,"StdHexExp::v_LaplacianMatrixOp_MatFree needs fixing");
/* if(mkey.GetNvariableLaplacianCoefficients() == 0)
{
// This implementation is only valid when there are no coefficients
// associated to the Laplacian operator
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nqtot = nquad0*nquad1;