Commit eb439808 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'feature/varopti_NEW' into 'master'

Add variational optimisation framework to NekMesh

This adds the variational mesh optimisation framework to NekMesh. Its now throughly tested for triangle, quad, tet and prism meshes along with triangle/quad and prism/tet hybrid meshes. Although not perfect we've got a pretty good grasp of how it behaves and why. Future work will be to move this code into a library as a "variational solver" where NekMesh utilises this along with alternative solvers.


See merge request !711
parents e025603b b3de3600
......@@ -76,6 +76,7 @@ v4.4.0
- Bug fix to get two meshgen regression tests working (!700)
- Remove libANN in deference to boost::geometry (!703)
- Refactor library to use NekMesh modules for CAD generation (!704)
- Add `varopti` process module to optimise meshes (!711)
- Add a mesh extract option to the linearise module to visualise the result
(!712)
- 2D to 3D mesh extrusion module (!715)
......
#SET(GraphSources
#SET(GraphSources
# GraphExample.cpp)
SET(MemoryManagerSources
SET(MemoryManagerSources
MemoryManager.cpp)
SET(PartitionAnalyseSources
PartitionAnalyse.cpp)
SET(FoundationSources
FoundationDemo.cpp)
FoundationDemo.cpp)
SET(NodalDemoSources NodalDemo.cpp)
SET(TimeIntegrationDemoSources
......
......@@ -151,10 +151,14 @@ int main(int argc, char *argv[])
{
util = new NodalUtilPrism(order, r, s, t);
}
else if(shape == eQuadrilateral)
{
util = new NodalUtilQuad(order, r, s);
}
ASSERTL1(util, "Unknown shape type!");
const int nPoints = r.num_elements();
const int dim = shape == eTriangle ? 2 : 3;
const int dim = (shape == eTriangle || shape == eQuadrilateral) ? 2 : 3;
if (vm.count("integral"))
{
......@@ -175,6 +179,9 @@ int main(int argc, char *argv[])
exact = -0.5 * (sin(1.0) + cos(1.0) + M_E * M_E *
(sin(1.0) - cos(1.0))) / M_E;
break;
case eQuadrilateral:
exact = 2.0 * (M_E - 1.0 / M_E) * sin(1.0);
break;
case eTetrahedron:
exact = 1.0 / M_E - 1.0 / M_E / M_E / M_E;
break;
......
......@@ -45,6 +45,7 @@ using namespace std;
#include <StdRegions/StdPrismExp.h>
#include <StdRegions/StdQuadExp.h>
#include <StdRegions/StdTetExp.h>
#include <StdRegions/StdHexExp.h>
#include <StdRegions/StdTriExp.h>
namespace Nektar
......@@ -295,6 +296,91 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
}
}
}
else if (geom->GetShapeType() == LibUtilities::eHexahedron)
{
vector<Array<OneD, NekDouble> > xyz;
for (int i = 0; i < geom->GetNumVerts(); i++)
{
Array<OneD, NekDouble> loc(3);
SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
p->GetCoords(loc);
xyz.push_back(loc);
}
Array<OneD, const LibUtilities::BasisSharedPtr> b = chi->GetBase();
Array<OneD, NekDouble> eta1 = b[0]->GetZ();
Array<OneD, NekDouble> eta2 = b[1]->GetZ();
Array<OneD, NekDouble> eta3 = b[2]->GetZ();
for (int k = 0; k < b[2]->GetNumPoints(); k++)
{
for (int j = 0; j < b[1]->GetNumPoints(); j++)
{
for (int i = 0; i < b[0]->GetNumPoints(); i++)
{
NekDouble a1 = 0.5 * (1 - eta1[i]);
NekDouble a2 = 0.5 * (1 + eta1[i]);
NekDouble b1 = 0.5 * (1 - eta2[j]),
b2 = 0.5 * (1 + eta2[j]);
NekDouble c1 = 0.5 * (1 - eta3[k]),
c2 = 0.5 * (1 + eta3[k]);
DNekMat dxdz(3, 3, 1.0, eFULL);
dxdz(0, 0) =
-0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0] +
0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0] -
0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0] +
0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
dxdz(1, 0) =
-0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1] +
0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1] -
0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1] +
0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
dxdz(2, 0) =
-0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2] +
0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2] -
0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2] +
0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
dxdz(0, 1) =
-0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0] +
0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0] -
0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0] +
0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
dxdz(1, 1) =
-0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1] +
0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1] -
0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1] +
0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
dxdz(2, 1) =
-0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2] +
0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2] -
0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2] +
0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
dxdz(0, 0) =
-0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0] -
0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0] +
0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0] +
0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
dxdz(1, 0) =
-0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1] -
0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1] +
0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1] +
0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
dxdz(2, 0) =
-0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2] -
0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2] +
0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2] +
0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
dxdz.Invert();
ret.push_back(dxdz);
}
}
}
}
else
{
ASSERTL0(false, "not coded");
......@@ -303,7 +389,9 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
return ret;
}
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(LocalRegions::ExpansionSharedPtr e, bool s)
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(
LocalRegions::ExpansionSharedPtr e,
bool s)
{
SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
StdRegions::StdExpansionSharedPtr chi = e->GetGeom()->GetXmap();
......@@ -359,6 +447,10 @@ Array<OneD, NekDouble> ProcessQualityMetric::GetQ(LocalRegions::ExpansionSharedP
chiMod = MemoryManager<StdRegions::StdPrismExp>::AllocateSharedPtr(
basisKeys[0], basisKeys[1], basisKeys[2]);
break;
case LibUtilities::eHexahedron:
chiMod = MemoryManager<StdRegions::StdHexExp>::AllocateSharedPtr(
basisKeys[0], basisKeys[1], basisKeys[2]);
break;
default:
ASSERTL0(false, "nope");
}
......
......@@ -140,6 +140,8 @@ SET(FoundationHeaders
./Foundations/NodalTetSPI.h
./Foundations/NodalPrismSPI.h
./Foundations/NodalTetSPIData.h
./Foundations/NodalQuadElec.h
./Foundations/NodalHexElec.h
./Foundations/NodalUtil.h
./Foundations/PhysGalerkinProject.h
./Foundations/Points.h
......@@ -166,6 +168,8 @@ SET(FoundationSources
./Foundations/NodalTriFekete.cpp
./Foundations/NodalTriSPI.cpp
./Foundations/NodalTetSPI.cpp
./Foundations/NodalQuadElec.cpp
./Foundations/NodalHexElec.cpp
./Foundations/NodalPrismSPI.cpp
./Foundations/NodalUtil.cpp
./Foundations/PhysGalerkinProject.cpp
......
......@@ -49,6 +49,8 @@
#include <LibUtilities/Foundations/NodalPrismSPI.h>
#include <LibUtilities/Foundations/NodalPrismEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalPrismElec.h>
#include <LibUtilities/Foundations/NodalQuadElec.h>
#include <LibUtilities/Foundations/NodalHexElec.h>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
......@@ -90,6 +92,8 @@ namespace Nektar
const bool NodalTriInited3 = PointsManager().RegisterCreator(PointsKey(0, eNodalTriSPI), NodalTriSPI::Create);
const bool NodalTriInited4 = PointsManager().RegisterCreator(PointsKey(0, eNodalTriEvenlySpaced), NodalTriEvenlySpaced::Create);
const bool NodalQuadInited1 = PointsManager().RegisterCreator(PointsKey(0, eNodalQuadElec), NodalQuadElec::Create);
const bool NodalTetInited1 = PointsManager().RegisterCreator(PointsKey(0, eNodalTetElec), NodalTetElec::Create);
const bool NodalTetInited2 = PointsManager().RegisterCreator(PointsKey(0, eNodalTetSPI), NodalTetSPI::Create);
const bool NodalTetInited3 = PointsManager().RegisterCreator(PointsKey(0, eNodalTetEvenlySpaced), NodalTetEvenlySpaced::Create);
......@@ -98,6 +102,8 @@ namespace Nektar
const bool NodalPrismInited2 = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismElec), NodalPrismElec::Create);
const bool NodalPrismInited3 = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismSPI), NodalPrismSPI::Create);
const bool NodalHexInited1 = PointsManager().RegisterCreator(PointsKey(0, eNodalHexElec), NodalHexElec::Create);
const bool Basis_Inited = BasisManager().RegisterGlobalCreator(Basis::Create);
};
......
///////////////////////////////////////////////////////////////////////////////
//
// File NodalHexElec.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: Nodal hexahedron with 3D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/NodalHexElec.h>
#include <LibUtilities/Foundations/NodalUtil.h>
#include <LibUtilities/Foundations/Points.h>
namespace Nektar
{
namespace LibUtilities
{
void NodalHexElec::CalculatePoints()
{
// Allocate the storage for points
unsigned int numPoints = GetNumPoints();
PointsKey e(numPoints, eGaussLobattoLegendre);
PointsManager()[e]->GetPoints(m_e0);
m_ew = PointsManager()[e]->GetW();
for (int i = 0; i < 3; i++)
{
m_points[i] = Array<OneD, DataType>(numPoints * numPoints * numPoints);
}
for (int k = 0, ct = 0; k < numPoints; k++)
{
for (int j = 0; j < numPoints; j++)
{
for (int i = 0; i < numPoints; i++, ct++)
{
m_points[0][ct] = m_e0[i];
m_points[1][ct] = m_e0[j];
m_points[2][ct] = m_e0[k];
}
}
}
}
void NodalHexElec::CalculateWeights()
{
unsigned int numPoints = GetNumPoints();
m_weights = Array<OneD, DataType>(numPoints * numPoints * numPoints);
for (int k = 0, ct = 0; k < numPoints; k++)
{
for (int j = 0; j < numPoints; j++)
{
for (int i = 0; i < numPoints; i++, ct++)
{
m_weights[ct] = m_ew[i] * m_ew[j] * m_ew[k];
}
}
}
}
void NodalHexElec::CalculateDerivMatrix()
{
}
boost::shared_ptr<PointsBaseType> NodalHexElec::Create(const PointsKey &key)
{
boost::shared_ptr<PointsBaseType> returnval(
MemoryManager<NodalHexElec>::AllocateSharedPtr(key));
returnval->Initialize();
return returnval;
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File NodalHexElec.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: Nodal hexahedron with 3D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NODALHEXELEC_H
#define NODALHEXELEC_H
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/Foundations/FoundationsFwd.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/LibUtilitiesDeclspec.h>
#include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
#include <boost/shared_ptr.hpp>
namespace Nektar
{
namespace LibUtilities
{
class NodalHexElec : public Points<NekDouble>
{
public:
virtual ~NodalHexElec()
{
}
LIB_UTILITIES_EXPORT static boost::shared_ptr<PointsBaseType> Create(
const PointsKey &key);
NodalHexElec(const PointsKey &key) : PointsBaseType(key)
{
}
private:
NodalHexElec() : PointsBaseType(NullPointsKey)
{
}
/// 1D GLL points.
Array<OneD, NekDouble> m_e0;
/// 1D GLL weights.
Array<OneD, NekDouble> m_ew;
void CalculatePoints();
void CalculateWeights();
void CalculateDerivMatrix();
};
} // end of namespace
} // end of namespace
#endif // NODALHEXELEC_H
///////////////////////////////////////////////////////////////////////////////
//
// File NodalQuadElec.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: Nodal quadrilateral with 2D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/NodalQuadElec.h>
#include <LibUtilities/Foundations/NodalUtil.h>
#include <LibUtilities/Foundations/Points.h>
namespace Nektar
{
namespace LibUtilities
{
void NodalQuadElec::CalculatePoints()
{
// Allocate the storage for points
unsigned int numPoints = GetNumPoints();
PointsKey e(numPoints, eGaussLobattoLegendre);
PointsManager()[e]->GetPoints(m_e0);
m_ew = PointsManager()[e]->GetW();
for (int i = 0; i < 2; i++)
{
m_points[i] = Array<OneD, DataType>(numPoints * numPoints);
}
for (int j = 0, ct = 0; j < numPoints; j++)
{
for (int i = 0; i < numPoints; i++, ct++)
{
m_points[0][ct] = m_e0[i];
m_points[1][ct] = m_e0[j];
}
}
}
void NodalQuadElec::CalculateWeights()
{
unsigned int numPoints = GetNumPoints();
m_weights = Array<OneD, DataType>(numPoints * numPoints);
for (int j = 0, ct = 0; j < numPoints; j++)
{
for (int i = 0; i < numPoints; i++, ct++)
{
m_weights[ct] = m_ew[i] * m_ew[j];
}
}
}
void NodalQuadElec::CalculateDerivMatrix()
{
}
boost::shared_ptr<PointsBaseType> NodalQuadElec::Create(const PointsKey &key)
{
boost::shared_ptr<PointsBaseType> returnval(
MemoryManager<NodalQuadElec>::AllocateSharedPtr(key));
returnval->Initialize();
return returnval;
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File NodalQuadElec.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: Nodal quadrilateral with 2D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NODALQUADELEC_H
#define NODALQUADELEC_H
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/Foundations/FoundationsFwd.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/LibUtilitiesDeclspec.h>
#include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
#include <boost/shared_ptr.hpp>
namespace Nektar
{
namespace LibUtilities
{
class NodalQuadElec : public Points<NekDouble>
{
public:
virtual ~NodalQuadElec()
{
}
LIB_UTILITIES_EXPORT static boost::shared_ptr<PointsBaseType> Create(
const PointsKey &key);
NodalQuadElec(const PointsKey &key) : PointsBaseType(key)
{
}
private:
NodalQuadElec() : PointsBaseType(NullPointsKey)
{
}
/// 1D GLL points
Array<OneD, NekDouble> m_e0;
/// 1D GLL weights
Array<OneD, NekDouble> m_ew;
void CalculatePoints();
void CalculateWeights();
void CalculateDerivMatrix();
};
} // end of namespace
} // end of namespace
#endif // NodalQuadElec
......@@ -967,11 +967,11 @@ NodalUtilHex::NodalUtilHex(int degree,
// Construct a mapping (i,j,k) -> m from the tensor product space (i,j,k) to
// a single ordering m.
for (int i = 0; i <= m_degree; ++i)
for (int k = 0; k <= m_degree; ++k)
{
for (int j = 0; j <= m_degree; ++j)
{
for (int k = 0; k <= m_degree - i; ++k)
for (int i = 0; i <= m_degree; ++i)
{
m_ordering.push_back(Mode(i, j, k));
}
......@@ -1053,11 +1053,18 @@ NekVector<NekDouble> NodalUtilHex::v_OrthoBasisDeriv(
ret[i] = jacobi_di[i] * jacobi_j[i] * jacobi_k[i];
}
}
else if (dir == 1)
{
for (int i = 0; i < m_numPoints; ++i)
{
ret[i] = jacobi_dj[i] * jacobi_i[i] * jacobi_k[i];
}
}
else
{
for (int i = 0; i < m_numPoints; ++i)
{
ret[i] = jacobi_i[i] * jacobi_dj[i] * jacobi_dk[i];
ret[i] = jacobi_i[i] * jacobi_j[i] * jacobi_dk[i];
}
}
......
......@@ -155,6 +155,8 @@ namespace Nektar
case eNodalTriElec:
case eNodalTriFekete:
case eNodalTriEvenlySpaced:
case eNodalTriSPI:
case eNodalQuadElec:
dimpoints = 2;
break;
......@@ -162,6 +164,7 @@ namespace Nektar
case eNodalTetEvenlySpaced:
case eNodalPrismEvenlySpaced:
case eNodalPrismElec:
case eNodalHexElec:
dimpoints = 3;
break;
......@@ -188,6 +191,10 @@ namespace Nektar
ASSERTL0(false,"this method cannot be implemented");
break;
case eNodalQuadElec:
totpoints = m_numpoints*m_numpoints;
break;
case eNodalTetElec:
case eNodalTetEvenlySpaced:
totpoints = m_numpoints*(m_numpoints+1)*(m_numpoints+2)/6;
......@@ -204,6 +211,10 @@ namespace Nektar
ASSERTL0(false,"this method cannot be implemented");
break;
case eNodalHexElec:
totpoints = m_numpoints*m_numpoints*m_numpoints;
break;
default:
break;
}
......
......@@ -68,7 +68,7 @@ namespace Nektar
eBoundaryLayerPoints, //!< 1D power law distribution for boundary layer points
eBoundaryLayerPointsRev, //!< 1D power law distribution for boundary layer points
eNodalTriElec, //!< 2D Nodal Electrostatic Points on a Triangle
eNodalTriFekete, //!< 2D Nodal Fekete Points on a Triangle
eNodalTriFekete, //!< 2D Nodal Fekete Points on a Triangle
eNodalTriEvenlySpaced, //!< 2D Evenly-spaced points on a Triangle
eNodalTetEvenlySpaced, //!< 3D Evenly-spaced points on a Tetrahedron
eNodalTetElec, //!< 3D Nodal Electrostatic Points on a Tetrahedron
......@@ -77,6 +77,8 @@ namespace Nektar
eNodalTriSPI, //!< 2D Nodal Symmetric positive internal triangle (Whitherden, Vincent)
eNodalTetSPI, //!< 3D Nodal Symmetric positive internal tet (Whitherden, Vincent)
eNodalPrismSPI, //!< 3D prism SPI
eNodalQuadElec, //!< 2D GLL for quad
eNodalHexElec, //!< 3D GLL for hex
SIZE_PointsType //