Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • nektar/redesign-prototypes
1 result
Show changes
Commits on Source (2)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1106 additions and 157 deletions
CMakeLists.txt
View file @ 229d7c18
......@@ -43,6 +43,7 @@ set(SRC Field.cpp Operators/Operator.cpp
Operators/OperatorDiagPrecon.cpp Operators/DiagPrecon/DiagPreconImpl.cpp
Operators/OperatorMatrix.cpp Operators/Matrix/MatrixImpl.cpp
Operators/OperatorNeuBndCond.cpp Operators/NeuBndCond/NeuBndCondImpl.cpp
Operators/OperatorRobBndCond.cpp Operators/RobBndCond/RobBndCondImpl.cpp
Operators/OperatorHelmholtz.cpp Operators/Helmholtz/HelmholtzImpl.cpp
Operators/OperatorIProductWRTDerivBase.cpp Operators/IProductWRTDerivBase/IProductWRTDerivBaseImpl.cpp
)
......
Operators/DirBndCond/DirBndCondStdMat.hpp
View file @ 229d7c18
......@@ -2,22 +2,14 @@
#include "Operators/OperatorDirBndCond.hpp"
#include <set>
#include <tuple>
#include <MultiRegions/AssemblyMap/AssemblyMapCG.h>
#include <MultiRegions/ContField.h>
using namespace Nektar;
using namespace Nektar::MultiRegions;
using namespace Nektar::SpatialDomains;
namespace Nektar::Operators::detail
{
void ImposeDirichletConditions(Array<OneD, NekDouble> &arr,
const ContFieldSharedPtr &contfield);
template <typename TData>
class OperatorDirBndCondImpl<TData, ImplStdMat>
: public OperatorDirBndCond<TData>
......@@ -28,106 +20,83 @@ public:
{
}
void apply(Field<TData, FieldState::Coeff> &inout) override
void apply(Field<TData, FieldState::Coeff> &out) override
{
// get number of local coeffs
auto contfield =
std::dynamic_pointer_cast<ContField>(this->m_expansionList);
auto nloc = contfield->GetLocalToGlobalMap()->GetNumLocalCoeffs();
// Field -> Array (** Needs changing!)
Array<OneD, NekDouble> arr(nloc, inout.GetStorage().GetCPUPtr());
// Core of function
ImposeDirichletConditions(arr, contfield);
// Array -> Field (** Needs changing!)
std::copy(arr.data(), arr.data() + nloc,
inout.GetStorage().GetCPUPtr());
}
// instantiation function for CreatorFunction in OperatorFactory
static std::unique_ptr<Operator<TData>> instantiate(
const MultiRegions::ExpListSharedPtr &expansionList)
{
return std::make_unique<OperatorDirBndCondImpl<TData, ImplStdMat>>(
expansionList);
}
// className - for OperatorFactory
static std::string className;
};
// Lifted code from Nektar (ContField.cpp)
void ImposeDirichletConditions(Array<OneD, NekDouble> &arr,
const ContFieldSharedPtr &contfield)
{
// get attributes from expansion list
auto locToGloMap = contfield->GetLocalToGlobalMap();
Array<OneD, ExpListSharedPtr> bndCondExpansions =
contfield->GetBndCondExpansions();
Array<OneD, BoundaryConditionShPtr> bndConditions =
contfield->GetBndConditions();
// lifted code:
int i, j;
int bndcnt = 0;
Array<OneD, NekDouble> sign =
locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
const Array<OneD, const int> map =
locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
for (i = 0; i < bndCondExpansions.size(); ++i)
{
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
auto &locToGloMap = contfield->GetLocalToGlobalMap();
auto &bndCondExpansions = contfield->GetBndCondExpansions();
auto &bndConditions = contfield->GetBndConditions();
auto &sign = locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
auto &map = locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
auto nloc = locToGloMap->GetNumLocalCoeffs();
size_t bndcnt = 0;
auto outptr = out.GetStorage().GetCPUPtr();
std::fill(outptr, outptr + nloc, 0.0);
for (size_t i = 0; i < bndCondExpansions.size(); ++i)
{
const Array<OneD, const NekDouble> bndcoeff =
(bndCondExpansions[i])->GetCoeffs();
if (locToGloMap->GetSignChange())
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
{
for (j = 0; j < (bndCondExpansions[i])->GetNcoeffs(); j++)
auto &bndcoeff = bndCondExpansions[i]->GetCoeffs();
if (locToGloMap->GetSignChange())
{
arr[map[bndcnt + j]] = sign[bndcnt + j] * bndcoeff[j];
for (size_t j = 0; j < bndCondExpansions[i]->GetNcoeffs();
j++)
{
*(outptr + map[bndcnt + j]) =
sign[bndcnt + j] * bndcoeff[j];
}
}
}
else
{
for (j = 0; j < (bndCondExpansions[i])->GetNcoeffs(); j++)
else
{
arr[map[bndcnt + j]] = bndcoeff[j];
for (size_t j = 0; j < bndCondExpansions[i]->GetNcoeffs();
j++)
{
*(outptr + map[bndcnt + j]) = bndcoeff[j];
}
}
}
bndcnt += bndCondExpansions[i]->GetNcoeffs();
}
bndcnt += bndCondExpansions[i]->GetNcoeffs();
}
// communicate local Dirichlet coeffs that are just
// touching a dirichlet boundary on another partition
std::set<int> &ParallelDirBndSign = locToGloMap->GetParallelDirBndSign();
// communicate local Dirichlet coeffs that are just
// touching a dirichlet boundary on another partition
auto &ParallelDirBndSign = locToGloMap->GetParallelDirBndSign();
for (auto &it : ParallelDirBndSign)
{
arr[it] *= -1;
}
for (auto &it : ParallelDirBndSign)
{
*(outptr + it) *= -1;
}
locToGloMap->UniversalAbsMaxBnd(arr);
Array<OneD, NekDouble> arr(nloc, outptr);
locToGloMap->UniversalAbsMaxBnd(arr);
std::copy(arr.get(), arr.get() + nloc, outptr);
for (auto &it : ParallelDirBndSign)
{
arr[it] *= -1;
for (auto &it : ParallelDirBndSign)
{
*(outptr + it) *= -1;
}
auto &copyLocalDirDofs = locToGloMap->GetCopyLocalDirDofs();
for (auto &it : copyLocalDirDofs)
{
*(outptr + std::get<0>(it)) =
*(outptr + std::get<1>(it)) * std::get<2>(it);
}
}
std::set<ExtraDirDof> &copyLocalDirDofs =
locToGloMap->GetCopyLocalDirDofs();
for (auto &it : copyLocalDirDofs)
// instantiation function for CreatorFunction in OperatorFactory
static std::unique_ptr<Operator<TData>> instantiate(
const MultiRegions::ExpListSharedPtr &expansionList)
{
arr[std::get<0>(it)] = arr[std::get<1>(it)] * std::get<2>(it);
return std::make_unique<OperatorDirBndCondImpl<TData, ImplStdMat>>(
expansionList);
}
}
// className - for OperatorFactory
static std::string className;
};
} // namespace Nektar::Operators::detail
Operators/NeuBndCond/NeuBndCondStdMat.hpp
View file @ 229d7c18
......@@ -10,10 +10,6 @@ using namespace Nektar::MultiRegions;
namespace Nektar::Operators::detail
{
void ImposeNeumannConditions(Array<OneD, NekDouble> &arr,
const ContFieldSharedPtr &contfield);
template <typename TData>
class OperatorNeuBndCondImpl<TData, ImplStdMat>
: public OperatorNeuBndCond<TData>
......@@ -26,20 +22,45 @@ public:
void apply(Field<TData, FieldState::Coeff> &inout) override
{
// get number of local coeffs
auto contField =
auto contfield =
std::dynamic_pointer_cast<ContField>(this->m_expansionList);
auto nloc = contField->GetLocalToGlobalMap()->GetNumLocalCoeffs();
// Field -> Array (** Needs changing!)
Array<OneD, NekDouble> rhs(nloc, inout.GetStorage().GetCPUPtr());
// Core of function
ImposeNeumannConditions(rhs, contField);
// Array -> Field (** Needs changing!)
std::copy(rhs.data(), rhs.data() + nloc,
inout.GetStorage().GetCPUPtr());
auto &locToGloMap = contfield->GetLocalToGlobalMap();
auto &bndCondExpansions = contfield->GetBndCondExpansions();
auto &bndConditions = contfield->GetBndConditions();
auto &sign = locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
auto &map = locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
size_t bndcnt = 0;
auto outptr = inout.GetStorage().GetCPUPtr();
// Add weak boundary conditions to forcing
for (size_t i = 0; i < bndCondExpansions.size(); ++i)
{
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eNeumann ||
bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eRobin)
{
auto &bndcoeff = bndCondExpansions[i]->GetCoeffs();
if (locToGloMap->GetSignChange())
{
for (size_t j = 0; j < bndCondExpansions[i]->GetNcoeffs();
j++)
{
*(outptr + map[bndcnt + j]) +=
sign[bndcnt + j] * bndcoeff[j];
}
}
else
{
for (size_t j = 0; j < bndCondExpansions[i]->GetNcoeffs();
j++)
{
*(outptr + map[bndcnt + j]) += bndcoeff[j];
}
}
}
bndcnt += bndCondExpansions[i]->GetNcoeffs();
}
}
// instantiation function for CreatorFunction in OperatorFactory
......@@ -54,50 +75,4 @@ public:
static std::string className;
};
// Lifted from ContField.cpp (HelmSolve function)
void ImposeNeumannConditions(Array<OneD, NekDouble> &arr,
const ContFieldSharedPtr &contfield)
{
auto locToGloMap = contfield->GetLocalToGlobalMap();
auto bndCondExpansions = contfield->GetBndCondExpansions();
auto bndConditions = contfield->GetBndConditions();
int bndcnt = 0;
Array<OneD, NekDouble> sign =
locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
const Array<OneD, const int> map =
locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
// Add weak boundary conditions to forcing
for (size_t i = 0; i < bndCondExpansions.size(); ++i)
{
if (bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eNeumann ||
bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eRobin)
{
const Array<OneD, const NekDouble> bndcoeff =
(bndCondExpansions[i])->GetCoeffs();
if (locToGloMap->GetSignChange())
{
for (size_t j = 0; j < (bndCondExpansions[i])->GetNcoeffs();
j++)
{
arr[map[bndcnt + j]] += sign[bndcnt + j] * bndcoeff[j];
}
}
else
{
for (size_t j = 0; j < (bndCondExpansions[i])->GetNcoeffs();
j++)
{
arr[map[bndcnt + j]] += bndcoeff[j];
}
}
}
bndcnt += bndCondExpansions[i]->GetNcoeffs();
}
}
} // namespace Nektar::Operators::detail
Operators/OperatorDirBndCond.hpp
View file @ 229d7c18
......@@ -19,7 +19,7 @@ public:
{
}
virtual void apply(Field<TData, FieldState::Coeff> &inout) = 0;
virtual void apply(Field<TData, FieldState::Coeff> &out) = 0;
};
// Descriptor / traits class for DirBndCond to be used by Operator create
......
Operators/OperatorRobBndCond.cpp 0 → 100644
View file @ 229d7c18
#include "OperatorRobBndCond.hpp"
namespace Nektar::Operators
{
// define static variables for RobBndCond operator descriptor for the default fp
// type config
template <> const std::string RobBndCond<default_fp_type>::key = "RobBndCond";
template <> const std::string RobBndCond<default_fp_type>::default_impl = "";
} // namespace Nektar::Operators
Operators/OperatorRobBndCond.hpp 0 → 100644
View file @ 229d7c18
#pragma once
#include "Field.hpp"
#include "Operator.hpp"
namespace Nektar::Operators
{
// Matrix operator base class
// Defines the apply operator to enforce apply parameter types
template <typename TData> class OperatorRobBndCond : public Operator<TData>
{
public:
virtual ~OperatorRobBndCond() = default;
OperatorRobBndCond(const MultiRegions::ExpListSharedPtr &expansionList)
: Operator<TData>(expansionList)
{
}
virtual void apply(Field<TData, FieldState::Coeff> &in,
Field<TData, FieldState::Coeff> &out,
const bool &negflag = false) = 0;
};
// Descriptor / traits class for RobBndCond to be used by Operator create
// function
template <typename TData = default_fp_type> struct RobBndCond
{
using class_name = OperatorRobBndCond<TData>;
static const std::string key;
static const std::string default_impl;
RobBndCond() = delete;
static std::shared_ptr<class_name> create(
const MultiRegions::ExpListSharedPtr &expansionList,
std::string pKey = "")
{
return Operator<TData>::template create<RobBndCond<TData>>(
expansionList, pKey);
}
};
namespace detail
{
// Template for implementation of RobBndCond operator
template <typename TData, typename Op> class OperatorRobBndCondImpl;
} // namespace detail
} // namespace Nektar::Operators
Operators/RobBndCond/RobBndCondImpl.cpp 0 → 100644
View file @ 229d7c18
#include "RobBndCondStdMat.hpp"
namespace Nektar::Operators::detail
{
// Register implementation with Operator Factory
template <>
std::string OperatorRobBndCondImpl<double, ImplStdMat>::className =
GetOperatorFactory<default_fp_type>().RegisterCreatorFunction(
"RobBndCond", OperatorRobBndCondImpl<double, ImplStdMat>::instantiate,
"");
} // namespace Nektar::Operators::detail
Operators/RobBndCond/RobBndCondStdMat.hpp 0 → 100644
View file @ 229d7c18
#pragma once
#include "Operators/OperatorRobBndCond.hpp"
#include <MultiRegions/AssemblyMap/AssemblyMapCG.h>
#include <MultiRegions/ContField.h>
using namespace Nektar;
using namespace Nektar::MultiRegions;
namespace Nektar::Operators::detail
{
template <typename TData>
class OperatorRobBndCondImpl<TData, ImplStdMat>
: public OperatorRobBndCond<TData>
{
public:
OperatorRobBndCondImpl(const MultiRegions::ExpListSharedPtr &expansionList)
: OperatorRobBndCond<TData>(std::move(expansionList))
{
}
void apply(Field<TData, FieldState::Coeff> &in,
Field<TData, FieldState::Coeff> &out,
const bool &negflag) override
{
auto ncoeffs = out.GetStorage().size();
auto *inptr = in.GetStorage().GetCPUPtr();
auto *outptr = out.GetStorage().GetCPUPtr();
Array<OneD, NekDouble> inarray(ncoeffs, inptr);
Array<OneD, NekDouble> robin(ncoeffs, 0.0);
auto *robptr = robin.get();
auto robinBCInfo = this->m_expansionList->GetRobinBCInfo();
for (auto &r : robinBCInfo) // add robin mass matrix
{
auto n = r.first;
auto offset = this->m_expansionList->GetCoeff_Offset(n);
auto expPtr = this->m_expansionList->GetExp(n);
Array<OneD, NekDouble> tmp;
for (auto rBC = r.second; rBC; rBC = rBC->next)
{
expPtr->AddRobinTraceContribution(
rBC->m_robinID, rBC->m_robinPrimitiveCoeffs,
inarray + offset, tmp = robin + offset);
}
}
if (negflag)
{
std::transform(robin.get(), robin.get() + ncoeffs, outptr, outptr,
[](const TData &rob, const TData &out)
{ return out - rob; });
}
else
{
std::transform(robin.get(), robin.get() + ncoeffs, outptr, outptr,
[](const TData &rob, const TData &out)
{ return out + rob; });
}
}
// instantiation function for CreatorFunction in OperatorFactory
static std::unique_ptr<Operator<TData>> instantiate(
const MultiRegions::ExpListSharedPtr &expansionList)
{
return std::make_unique<OperatorRobBndCondImpl<TData, ImplStdMat>>(
expansionList);
}
// className - for OperatorFactory
static std::string className;
};
} // namespace Nektar::Operators::detail
run/Helmholtz2D_P7_AllBCs.xml 0 → 100644
View file @ 229d7c18
<?xml version="1.0" encoding="utf-8"?>
<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
<GEOMETRY DIM="2" SPACE="2">
<VERTEX>
<V ID="0"> -1.000000000000000 3.500000000000000 0.0 </V>
<V ID="1"> -1.000000000000000 0.500000000000000 0.0 </V>
<V ID="2"> -1.000000000000000 2.500000000000000 0.0 </V>
<V ID="3"> -1.000000000000000 1.500000000000000 0.0 </V>
<V ID="4"> 3.800000000000000 4.500000000000000 0.0 </V>
<V ID="5"> 0.200000000000000 4.500000000000000 0.0 </V>
<V ID="6"> 2.900000000000000 4.500000000000000 0.0 </V>
<V ID="7"> 2.000000000000000 4.500000000000000 0.0 </V>
<V ID="8"> 1.100000000000000 4.500000000000000 0.0 </V>
<V ID="9"> 5.000000000000000 0.500000000000000 0.0 </V>
<V ID="10"> 5.000000000000000 3.500000000000000 0.0 </V>
<V ID="11"> 5.000000000000000 1.500000000000000 0.0 </V>
<V ID="12"> 5.000000000000000 2.500000000000000 0.0 </V>
<V ID="13"> 0.200000000000000 -0.500000000000000 0.0 </V>
<V ID="14"> 3.800000000000000 -0.500000000000000 0.0 </V>
<V ID="15"> 1.100000000000000 -0.500000000000000 0.0 </V>
<V ID="16"> 2.000000000000000 -0.500000000000000 0.0 </V>
<V ID="17"> 2.900000000000000 -0.500000000000000 0.0 </V>
<V ID="18"> -0.400000000000000 4.000000000000000 0.0 </V>
<V ID="19"> 4.400000000000000 4.000000000000000 0.0 </V>
<V ID="20"> 4.400000000000000 0.0 0.0 </V>
<V ID="21"> -0.400000000000000 0.0 0.0 </V>
<V ID="22"> -0.040000000000000 2.700000000000000 0.0 </V>
<V ID="23"> 0.920000000000000 1.900000000000000 0.0 </V>
<V ID="24"> 1.880000000000000 1.100000000000000 0.0 </V>
<V ID="25"> 2.840000000000000 0.300000000000000 0.0 </V>
<V ID="26"> -0.119314370713000 1.785562178050000 0.0 </V>
<V ID="27"> 1.713659159880000 0.169773145192000 0.0 </V>
<V ID="28"> 0.677713625957000 1.180143861910000 0.0 </V>
<V ID="29"> -0.188491169544000 0.869785275722000 0.0 </V>
<V ID="30"> 0.715978655871000 0.336366878402000 0.0 </V>
<V ID="31"> 3.289084972900000 0.882472796343000 0.0 </V>
<V ID="32"> 3.682080702820000 1.501561993610000 0.0 </V>
<V ID="33"> 3.972007805980000 2.179240079990000 0.0 </V>
<V ID="34"> 4.218223358320000 2.817871076440000 0.0 </V>
<V ID="35"> 4.379282019760000 3.315489806910000 0.0 </V>
<V ID="36"> 3.668892462910000 3.970508359810000 0.0 </V>
<V ID="37"> 3.155413025260000 3.802413365130000 0.0 </V>
<V ID="38"> 2.269709494930000 3.615953163480000 0.0 </V>
<V ID="39"> 1.341308209110000 3.906056835480000 0.0 </V>
<V ID="40"> 0.934565531168000 3.565834210030000 0.0 </V>
<V ID="41"> 0.461690231830000 3.150972034220000 0.0 </V>
<V ID="42"> 1.383772742880000 2.431995597660000 0.0 </V>
<V ID="43"> 2.324656705230000 1.661732140030000 0.0 </V>
<V ID="44"> 3.642425028600000 3.246403386770000 0.0 </V>
<V ID="45"> 4.022057777820000 3.634175940310000 0.0 </V>
<V ID="46"> 1.833157069900000 2.985307743670000 0.0 </V>
<V ID="47"> 2.755342831120000 2.226765404480000 0.0 </V>
<V ID="48"> 3.179196984040000 2.779077508740000 0.0 </V>
</VERTEX>
<EDGE>
<E ID="0"> 1 21 </E>
<E ID="1"> 13 21 </E>
<E ID="2"> 13 15 </E>
<E ID="3"> 15 16 </E>
<E ID="4"> 16 17 </E>
<E ID="5"> 14 17 </E>
<E ID="6"> 1 3 </E>
<E ID="7"> 1 29 </E>
<E ID="8"> 21 29 </E>
<E ID="9"> 21 30 </E>
<E ID="10"> 13 30 </E>
<E ID="11"> 15 30 </E>
<E ID="12"> 15 27 </E>
<E ID="13"> 16 27 </E>
<E ID="14"> 16 25 </E>
<E ID="15"> 17 25 </E>
<E ID="16"> 3 29 </E>
<E ID="17"> 29 30 </E>
<E ID="18"> 27 30 </E>
<E ID="19"> 25 27 </E>
<E ID="20"> 2 3 </E>
<E ID="21"> 3 26 </E>
<E ID="22"> 26 29 </E>
<E ID="23"> 28 29 </E>
<E ID="24"> 28 30 </E>
<E ID="25"> 24 30 </E>
<E ID="26"> 24 27 </E>
<E ID="27"> 2 26 </E>
<E ID="28"> 26 28 </E>
<E ID="29"> 24 28 </E>
<E ID="30"> 0 2 </E>
<E ID="31"> 2 22 </E>
<E ID="32"> 22 26 </E>
<E ID="33"> 23 26 </E>
<E ID="34"> 23 28 </E>
<E ID="35"> 0 22 </E>
<E ID="36"> 22 23 </E>
<E ID="37"> 23 24 </E>
<E ID="38"> 24 25 </E>
<E ID="39"> 14 25 </E>
<E ID="40"> 0 18 </E>
<E ID="41"> 22 41 </E>
<E ID="42"> 23 42 </E>
<E ID="43"> 24 43 </E>
<E ID="44"> 25 31 </E>
<E ID="45"> 14 20 </E>
<E ID="46"> 18 41 </E>
<E ID="47"> 41 42 </E>
<E ID="48"> 42 43 </E>
<E ID="49"> 31 43 </E>
<E ID="50"> 20 31 </E>
<E ID="51"> 5 18 </E>
<E ID="52"> 40 41 </E>
<E ID="53"> 42 46 </E>
<E ID="54"> 43 47 </E>
<E ID="55"> 31 32 </E>
<E ID="56"> 9 20 </E>
<E ID="57"> 5 40 </E>
<E ID="58"> 40 46 </E>
<E ID="59"> 46 47 </E>
<E ID="60"> 32 47 </E>
<E ID="61"> 9 32 </E>
<E ID="62"> 5 8 </E>
<E ID="63"> 39 40 </E>
<E ID="64"> 38 46 </E>
<E ID="65"> 47 48 </E>
<E ID="66"> 32 33 </E>
<E ID="67"> 9 11 </E>
<E ID="68"> 8 39 </E>
<E ID="69"> 7 8 </E>
<E ID="70"> 38 39 </E>
<E ID="71"> 7 38 </E>
<E ID="72"> 38 48 </E>
<E ID="73"> 33 48 </E>
<E ID="74"> 11 33 </E>
<E ID="75"> 6 7 </E>
<E ID="76"> 37 38 </E>
<E ID="77"> 44 48 </E>
<E ID="78"> 33 34 </E>
<E ID="79"> 11 12 </E>
<E ID="80"> 6 37 </E>
<E ID="81"> 37 44 </E>
<E ID="82"> 34 44 </E>
<E ID="83"> 12 34 </E>
<E ID="84"> 4 6 </E>
<E ID="85"> 36 37 </E>
<E ID="86"> 44 45 </E>
<E ID="87"> 34 35 </E>
<E ID="88"> 10 12 </E>
<E ID="89"> 4 36 </E>
<E ID="90"> 36 45 </E>
<E ID="91"> 35 45 </E>
<E ID="92"> 10 35 </E>
<E ID="93"> 19 45 </E>
<E ID="94"> 4 19 </E>
<E ID="95"> 10 19 </E>
</EDGE>
<ELEMENT>
<T ID="0"> 6 7 16 </T>
<T ID="1"> 0 8 7 </T>
<T ID="2"> 9 17 8 </T>
<T ID="3"> 1 10 9 </T>
<T ID="4"> 2 11 10 </T>
<T ID="5"> 11 12 18 </T>
<T ID="6"> 3 13 12 </T>
<T ID="7"> 14 19 13 </T>
<T ID="8"> 4 15 14 </T>
<T ID="9"> 5 39 15 </T>
<T ID="10"> 21 27 20 </T>
<T ID="11"> 16 22 21 </T>
<T ID="12"> 23 28 22 </T>
<T ID="13"> 17 24 23 </T>
<T ID="14"> 24 25 29 </T>
<T ID="15"> 18 26 25 </T>
<T ID="16"> 19 38 26 </T>
<T ID="17"> 30 31 35 </T>
<T ID="18"> 27 32 31 </T>
<T ID="19"> 32 33 36 </T>
<T ID="20"> 28 34 33 </T>
<T ID="21"> 29 37 34 </T>
<Q ID="22"> 35 41 46 40 </Q>
<Q ID="23"> 36 42 47 41 </Q>
<Q ID="24"> 37 43 48 42 </Q>
<Q ID="25"> 38 44 49 43 </Q>
<Q ID="26"> 39 45 50 44 </Q>
<Q ID="27"> 46 52 57 51 </Q>
<Q ID="28"> 47 53 58 52 </Q>
<Q ID="29"> 48 54 59 53 </Q>
<Q ID="30"> 49 55 60 54 </Q>
<Q ID="31"> 50 56 61 55 </Q>
<Q ID="32"> 57 63 68 62 </Q>
<Q ID="33"> 58 64 70 63 </Q>
<Q ID="34"> 59 65 72 64 </Q>
<Q ID="35"> 60 66 73 65 </Q>
<Q ID="36"> 61 67 74 66 </Q>
<Q ID="37"> 68 70 71 69 </Q>
<Q ID="38"> 71 76 80 75 </Q>
<Q ID="39"> 72 77 81 76 </Q>
<Q ID="40"> 73 78 82 77 </Q>
<Q ID="41"> 74 79 83 78 </Q>
<Q ID="42"> 80 85 89 84 </Q>
<Q ID="43"> 81 86 90 85 </Q>
<Q ID="44"> 82 87 91 86 </Q>
<Q ID="45"> 83 88 92 87 </Q>
<Q ID="46"> 89 90 93 94 </Q>
<Q ID="47"> 91 92 95 93 </Q>
</ELEMENT>
<COMPOSITE>
<C ID="0"> Q[22-47] </C>
<C ID="1"> T[0-21] </C>
<C ID="2"> E[0-1] </C>
<C ID="3"> E[2-5] </C>
<C ID="4"> E[45] </C>
<C ID="5"> E[56] </C>
<C ID="6"> E[67] </C>
<C ID="7"> E[79] </C>
<C ID="8"> E[88] </C>
<C ID="9"> E[94-95] </C>
<C ID="10"> E[84,75,69,62,51,40,30,20,6] </C>
</COMPOSITE>
<DOMAIN> C[0-1] </DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
<E COMPOSITE="C[1]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<SOLVERINFO>
<I PROPERTY="GLOBALSYSSOLN" VALUE="IterativeFull" />
<I PROPERTY="LinSysIterSolver" VALUE="ConjugateGradientLoc"/>
<I PROPERTY="Preconditioner" VALUE="Diagonal"/>
</SOLVERINFO>
<PARAMETERS>
<P> Lambda = 1 </P>
<P> IterativeSolverTolerance = 1e-14 </P>
</PARAMETERS>
<VARIABLES>
<V ID="0"> u </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[2] </B>
<B ID="1"> C[3] </B>
<B ID="2"> C[4] </B>
<B ID="3"> C[5] </B>
<B ID="4"> C[6] </B>
<B ID="5"> C[7] </B>
<B ID="6"> C[8] </B>
<B ID="7"> C[9] </B>
<B ID="8"> C[10] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</REGION>
<REGION REF="1">
<R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1-PI*sin(PI*x)*cos(PI*y)"
PRIMCOEFF="1" />
</REGION>
<REGION REF="2">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</REGION>
<REGION REF="3">
<N VAR="u"
VALUE="(5/sqrt(61))*PI*cos(PI*x)*sin(PI*y)-(6/sqrt(61))*PI*sin(PI*x)*cos(PI*y)" />
</REGION>
<REGION REF="4">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</REGION>
<REGION REF="5">
<N VAR="u" VALUE="PI*cos(PI*x)*sin(PI*y)" />
</REGION>
<REGION REF="6">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</REGION>
<REGION REF="7">
<R VAR="u"
VALUE="sin(PI*x)*sin(PI*y) -1+ (5/sqrt(61))*PI*cos(PI*x)*sin(PI*y)+(6/sqrt(61))*PI*sin(PI*x)*cos(PI*y)"
PRIMCOEFF="1" />
</REGION>
<REGION REF="8">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Forcing">
<E VAR="u" VALUE="-Lambda*(sin(PI*x)*sin(PI*y)-1) - 2*PI*PI*sin(PI*x)*sin(PI*y)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="sin(PI*x)*sin(PI*y)-1" />
</FUNCTION>
</CONDITIONS>
</NEKTAR>
run/Helmholtz3D_Hex_AllBCs_P6.xml 0 → 100644
View file @ 229d7c18
This diff is collapsed.
tests/CMakeLists.txt
View file @ 229d7c18
......@@ -185,6 +185,15 @@ target_include_directories(test_assmbscatr PRIVATE ${NEKTAR++_INCLUDE_DIRS})
target_compile_definitions(test_assmbscatr PRIVATE
-DTEST_PATH="${CMAKE_SOURCE_DIR}")
add_executable(test_dirichlet test_dirichlet.cpp)
target_link_libraries(test_dirichlet PRIVATE Operators)
target_link_libraries(test_dirichlet PRIVATE ${NEKTAR++_LIBRARIES})
target_link_libraries(test_dirichlet PRIVATE Boost::unit_test_framework)
target_include_directories(test_dirichlet PRIVATE "${CMAKE_SOURCE_DIR}")
target_include_directories(test_dirichlet PRIVATE ${NEKTAR++_INCLUDE_DIRS})
target_compile_definitions(test_dirichlet PRIVATE
-DTEST_PATH="${CMAKE_SOURCE_DIR}")
add_executable(test_fwdtrans test_fwdtrans.cpp)
target_link_libraries(test_fwdtrans PRIVATE Operators)
target_link_libraries(test_fwdtrans PRIVATE ${NEKTAR++_LIBRARIES})
......@@ -340,6 +349,12 @@ add_test(
WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
)
add_test(
NAME Dirichlet
COMMAND test_dirichlet
WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
)
add_test(
NAME FwdTrans
COMMAND test_fwdtrans
......
tests/init_dirichletfields.hpp 0 → 100644
View file @ 229d7c18
#include "init_fields.hpp"
using namespace std;
using namespace Nektar::Operators;
using namespace Nektar::LibUtilities;
using namespace Nektar;
class DirichletField
: public InitFields<double, FieldState::Coeff, FieldState::Coeff,
MultiRegions::ContField>
{
public:
DirichletField()
: InitFields<double, FieldState::Coeff, FieldState::Coeff,
MultiRegions::ContField>()
{
}
void ExpectedSolution(const std::vector<BlockAttributes> &blocks,
double *inptr)
{
Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
// Calculate expected result from Nektar++
fixt_explist->ImposeDirichletConditions(outcoeffs);
// Copy expected result from Array to fixt_expected
double *coeffptr = outcoeffs.get();
for (auto const &block : blocks)
{
for (size_t el = 0; el < block.num_elements; ++el)
{
for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
{
(*inptr++) = (*coeffptr++);
}
}
for (size_t el = 0; el < block.num_padding_elements; ++el)
{
for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
{
inptr++;
}
}
}
}
};
class Helmholtz2D_P7_AllBCs : public DirichletField
{
public:
Helmholtz2D_P7_AllBCs()
{
meshName = "run/Helmholtz2D_P7_AllBCs.xml";
}
};
class Helmholtz3D_Hex_AllBCs_P6 : public DirichletField
{
public:
Helmholtz3D_Hex_AllBCs_P6()
{
meshName = "run/Helmholtz3D_Hex_AllBCs_P6.xml";
}
};
tests/init_fields.hpp
View file @ 229d7c18
......@@ -92,7 +92,7 @@ public:
{
fixt_explist =
MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
session, graph, "DefaultVar", true, false,
session, graph, "u", true, false,
Collections::eNoCollection);
}
......@@ -100,8 +100,7 @@ public:
{
fixt_explist =
MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
session, graph, true, "DefaultVar",
Collections::eNoCollection);
session, graph, true, "u", Collections::eNoCollection);
}
// Generate a blocks definition from the expansion list for each state
......
tests/test_dirichlet.cpp 0 → 100644
View file @ 229d7c18
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE TestDirichlet
#include <boost/test/tools/output_test_stream.hpp>
#include <boost/test/unit_test.hpp>
#include <iostream>
#include <memory>
#include "Operators/OperatorDirBndCond.hpp"
#include "init_dirichletfields.hpp"
BOOST_AUTO_TEST_SUITE(TestDirichlet)
using namespace std;
using namespace Nektar::Operators;
using namespace Nektar::LibUtilities;
using namespace Nektar;
BOOST_FIXTURE_TEST_CASE(dirichlet2d_p7_allbcs, Helmholtz2D_P7_AllBCs)
{
Configure();
auto DirBndCondOp = DirBndCond<>::create(fixt_explist, "");
DirBndCondOp->apply(*fixt_out);
ExpectedSolution(fixt_expected->GetBlocks(),
fixt_expected->GetStorage().GetCPUPtr());
BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
boost::test_tools::output_test_stream output;
{
OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
}
}
BOOST_FIXTURE_TEST_CASE(dirichlet3d_hex_allbcs_p6, Helmholtz3D_Hex_AllBCs_P6)
{
Configure();
auto DirBndCondOp = DirBndCond<>::create(fixt_explist, "");
DirBndCondOp->apply(*fixt_out);
ExpectedSolution(fixt_expected->GetBlocks(),
fixt_expected->GetStorage().GetCPUPtr());
BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
boost::test_tools::output_test_stream output;
{
OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
}
}
BOOST_AUTO_TEST_SUITE_END()