Select Git revision
LocalDtAcousticSolver.cpp
LocalDtAcousticSolver.cpp 47.85 KiB
#include <scheme/LocalDtAcousticSolver.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/ItemValueVariant.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
#include <variant>
#include <vector>
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
: public LocalDtAcousticSolverHandler::ILocalDtAcousticSolver
{
private:
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
using MeshType = Mesh<Connectivity<Dimension>>;
using MeshDataType = MeshData<Dimension>;
using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
class FixedBoundaryCondition;
class PressureBoundaryCondition;
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
class ExternalFSIVelocityBoundaryCondition;
class CouplingBoundaryCondition;
using BoundaryCondition = std::variant<FixedBoundaryCondition,
PressureBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition,
ExternalFSIVelocityBoundaryCondition,
CouplingBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
NodeValuePerCell<const Rdxd>
_computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValuePerCell<const Rd> njr = mesh_data.njr();
NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const size_t& nb_nodes = Ajr.numberOfSubValues(j);
const double& rhoc_j = rhoc[j];
for (size_t r = 0; r < nb_nodes; ++r) {
Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r));
} });
return Ajr;
}
NodeValuePerCell<const Rdxd>
_computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
parallel_for(
Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
const auto& cell_faces = cell_to_face_matrix[j];
const double& rho_c = rhoc[j];
for (size_t L = 0; L < cell_faces.size(); ++L) {
const FaceId& l = cell_faces[L];
const auto& face_nodes = face_to_node_matrix[l];
auto local_node_number_in_cell = [&](NodeId node_number) {
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
if (node_number == cell_nodes[i_node]) {
return i_node;
}
}
return std::numeric_limits<size_t>::max();
};
for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
const size_t R = local_node_number_in_cell(face_nodes[rl]);
Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
}
}
});
return Ajr;
}
NodeValuePerCell<const Rdxd>
_computeAjr(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
if constexpr (Dimension == 1) {
return _computeGlaceAjr(mesh, rhoc);
} else {
switch (solver_type) {
case SolverType::Glace: {
return _computeGlaceAjr(mesh, rhoc);
}
case SolverType::Eucclhyd: {
return _computeEucclhydAjr(mesh, rhoc);
}
default: {
throw UnexpectedError("invalid solver type");
}
} }
}
NodeValue<Rdxd>
_computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
{
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rdxd> Ar{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
Rdxd sum = zero;
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
for (size_t j = 0; j < node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
sum += Ajr(J, R);
}
Ar[r] = sum;
});
return Ar;
}
NodeValue<Rd>
_computeBr(const Mesh<Connectivity<Dimension>>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rd> b{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
Rd br = zero;
for (size_t j = 0; j < node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R);
}
b[r] = br;
});
return b;
}
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
BoundaryConditionList bc_list;
for (const auto& bc_descriptor : bc_descriptor_list) {
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::fixed: {
bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::external: {
if constexpr (Dimension == 1) {
const ExternalBoundaryConditionDescriptor& extern_bc_descriptor =
dynamic_cast<const ExternalBoundaryConditionDescriptor&>(*bc_descriptor);
bc_list.emplace_back(
ExternalFSIVelocityBoundaryCondition(mesh, getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
extern_bc_descriptor.socket()));
} else {
throw NotImplementedError("external FSI velocity not implemented for dimension > 1");
}
break;
}
case IBoundaryConditionDescriptor::Type::dirichlet: {
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if (dirichlet_bc_descriptor.name() == "velocity") {
MeshNodeBoundary<Dimension> mesh_node_boundary =
getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
Array<const Rd> value_list =
InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
mesh.xr(),
mesh_node_boundary.nodeList());
bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
} else if (dirichlet_bc_descriptor.name() == "pressure") {
const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
if constexpr (Dimension == 1) {
MeshNodeBoundary<Dimension> mesh_node_boundary =
getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
Array<const double> node_values =
InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
mesh_node_boundary.nodeList());
bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
} else {
MeshFaceBoundary<Dimension> mesh_face_boundary =
getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
Array<const double> face_values =
InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
mesh_face_boundary.faceList());
bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
}
} else {
is_valid_boundary_condition = false;
}
break;
}
case IBoundaryConditionDescriptor::Type::coupling: {
const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
bc_list.emplace_back(
CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),coupling_bc_descriptor.label()));
break;
} default: {
is_valid_boundary_condition = false;
}
}
if (not is_valid_boundary_condition) {
std::ostringstream error_msg;
error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
throw NormalError(error_msg.str());
}
}
return bc_list;
}
void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyExternalVelocityBC(const BoundaryConditionList& bc_list,
const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const;
void _applyCouplingBC(NodeValue<Rdxd>& Ar1,
NodeValue<Rdxd>& Ar2,
NodeValue<Rd>& br1,
NodeValue<Rd>& br2,
const std::vector<NodeId>& map1,
const std::vector<NodeId>& map2) const;
void _applyCouplingBC(const MeshType& mesh,
NodeValue<Rd>& ur,
NodeValue<Rd>& CR_ur,
NodeValuePerCell<Rd>& Fjr,
NodeValuePerCell<Rd>& CR_Fjr,
const std::vector<NodeId>& map2) const;
void
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
this->_applyPressureBC(bc_list, mesh, br);
this->_applySymmetryBC(bc_list, Ar, br);
this->_applyVelocityBC(bc_list, Ar, br);
}
NodeValue<Rd>
_computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
NodeValue<Rd> u{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
return u;
}
NodeValuePerCell<Rd>
_computeFjr(const MeshType& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const NodeValue<const Rd>& ur,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValuePerCell<Rd> F{mesh.connectivity()};
parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
for (size_t r = 0; r < cell_nodes.size(); ++r) {
F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
}
});
return F;
}
std::tuple<std::vector<NodeId>,
std::vector<NodeId>>
_computeMapping(std::shared_ptr<const IMesh>& i_mesh1,
std::shared_ptr<const IMesh>& i_mesh2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
std::vector<NodeId> map1;
std::vector<NodeId> map2;
NodeValue<Rd> xr1 = copy(mesh1.xr());
NodeValue<Rd> xr2 = copy(mesh2.xr());
for(const auto& boundary_condition1 : bc_list1){
std::visit([&](auto&& bc1) {
using T1 = std::decay_t<decltype(bc1)>;
if constexpr (std::is_same_v<CouplingBoundaryCondition,T1>){
const auto& node_list1 = bc1.nodeList();
for(const auto& boundary_condition2 : bc_list2){
std::visit([&](auto&& bc2) {
using T2 = std::decay_t<decltype(bc2)>;
if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
if(bc1.label() == bc2.label()){
const auto& node_list2 = bc2.nodeList();
for(size_t i = 0; i<node_list1.size(); i++){
NodeId node_id1 = node_list1[i];
NodeId node_id2;
for(size_t j = 0; j<node_list2.size(); j++){
node_id2 = node_list2[j];
double err = 0;
for(size_t j = 0; j<Dimension; j++){
err += (xr1[node_id1][j] - xr2[node_id2][j])*(xr1[node_id1][j] - xr2[node_id2][j]);
}
if(sqrt(err) < 1e-14){
map1.push_back(node_id1);
map2.push_back(node_id2);
}
}
}
}
}
},boundary_condition2);
}
}
},boundary_condition1);
}
return std::make_tuple(map1,map2);
}
public:
std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v, const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
if (not i_mesh) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
this->_applyExternalVelocityBC(bc_list, p, Ar, br);
synchronize(Ar);
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
std::tuple<const std::shared_ptr<const ItemValueVariant>,
const std::shared_ptr<const SubItemValuePerItemVariant>,
const std::shared_ptr<const ItemValueVariant>,
const std::shared_ptr<const SubItemValuePerItemVariant>,
NodeValue<Rd>,
NodeValuePerCell<Rd>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
const std::vector<NodeId>& map1,
const std::vector<NodeId>& map2) const
{
std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
if (not i_mesh1) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho1_v, c1_v, u1_v, p1_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
if (not i_mesh2) {
throw NormalError("discrete functions are not defined on the same mesh"); }
if (not checkDiscretizationType({rho2_v, c2_v, u2_v, p2_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c1 = c1_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p1 = p1_v->get<DiscreteScalarFunction>();
const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c2 = c2_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p2 = p2_v->get<DiscreteScalarFunction>();
NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * c1);
NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * c2);
NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, p1);
NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, p2);
const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
this->_applyCouplingBC(Ar1,Ar2,br1,br2,map1,map2);
synchronize(Ar1);
synchronize(br1);
synchronize(Ar2);
synchronize(br2);
NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
std::make_shared<const ItemValueVariant>(ur2),
std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
CR_ur,
CR_Fjr);
}
std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
NodeValue<Rd> CR_ur,
NodeValuePerCell<Rd> CR_Fjr,
const std::vector<NodeId> map2) const
{
std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
if (not i_mesh) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
this->_applyExternalVelocityBC(bc_list, p, Ar, br);
synchronize(Ar);
synchronize(br);
NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const NodeValue<const Rd>& ur,
const NodeValuePerCell<const Rd>& Fjr) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) { const NodeId r = cell_nodes[R];
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
new_u[j] -= dt_over_Mj * momentum_fluxes;
new_E[j] -= dt_over_Mj * energy_fluxes;
});
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
}
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply_fluxes(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
{
std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
if (not i_mesh) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
return this->apply_fluxes(dt, //
dynamic_cast<const MeshType&>(*i_mesh), //
rho_v->get<DiscreteScalarFunction>(), //
u_v->get<DiscreteVectorFunction>(), //
E_v->get<DiscreteScalarFunction>(), //
ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
mesh_correction(const MeshType& mesh1,
const MeshType& mesh2,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
std::vector<NodeId>& map1,
std::vector<NodeId>& map2) const
{
NodeValue<Rd> xr1 = copy(mesh1.xr());
NodeValue<Rd> new_xr2 = copy(mesh2.xr());
size_t n = map1.size();
for(size_t i = 0; i<n; i++){
for(size_t j = 0; j<Dimension; j++){
new_xr2[map2[i]][j] = xr1[map1[i]][j];
} }
std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
return {new_mesh2,std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh2, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_E))};
}
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
mesh_correction(const std::shared_ptr<const IMesh>& i_mesh1,
const std::shared_ptr<const IMesh>& i_mesh2,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
std::vector<NodeId>& map1,
std::vector<NodeId>& map2) const
{
return this->mesh_correction(dynamic_cast<const MeshType&>(*i_mesh1),
dynamic_cast<const MeshType&>(*i_mesh2),
rho->get<DiscreteScalarFunction>(),
u->get<DiscreteVectorFunction>(),
E->get<DiscreteScalarFunction>(),
map1,
map2);
}
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
const double& dt1,
const size_t& q,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
const std::shared_ptr<const DiscreteFunctionVariant>& u1,
const std::shared_ptr<const DiscreteFunctionVariant>& u2,
const std::shared_ptr<const DiscreteFunctionVariant>& E1,
const std::shared_ptr<const DiscreteFunctionVariant>& E2,
const std::shared_ptr<const DiscreteFunctionVariant>& c1,
const std::shared_ptr<const DiscreteFunctionVariant>& c2,
const std::shared_ptr<const DiscreteFunctionVariant>& p1,
const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
std::shared_ptr<const IMesh> new_m2 = getCommonMesh({rho2, c2, u2, p2});
std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2,map1,map2);
const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
double dt2 = dt1/q;
const double& gamma = 1.4;
double sum_dt = 0;
do{
const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 = std::make_shared<const DiscreteFunctionVariant>(p_d);
const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 = std::make_shared<const DiscreteFunctionVariant>(c_d);
auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
sum_dt += dt2;
}while(sum_dt < dt1);
std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
}
LocalDtAcousticSolver() = default;
LocalDtAcousticSolver(LocalDtAcousticSolver&&) = default;
~LocalDtAcousticSolver() = default;
};
template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
Assert(node_cell_list.size() == 1);
CellId node_cell_id = node_cell_list[0];
size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
});
} else {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
NodeId node_id = face_nodes[i_node];
br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
}
}
}
}
},
boundary_condition);
}
}
template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
const Rd& n = bc.outgoingNormal();
const Rdxd I = identity;
const Rdxd nxn = tensorProduct(n, n);
const Rdxd P = I - nxn;
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
const NodeId r = node_list[r_number];
Ar[r] = P * Ar[r] * P + nxn;
br[r] = P * br[r];
});
}
},
boundary_condition);
}
}
template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
const auto& value = value_list[i_node];
Ar[node_id] = identity;
br[node_id] = value;
});
} else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
Ar[node_id] = identity;
br[node_id] = zero;
});
}
},
boundary_condition);
}
}
template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVelocityBC(
const BoundaryConditionList& bc_list,
const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<ExternalFSIVelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList(p);
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
const auto& value = value_list[i_node];
Ar[node_id] = identity;
br[node_id] = value;
});
}
},
boundary_condition);
}
}
template <size_t Dimension>
void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(NodeValue<Rdxd>& Ar1,
NodeValue<Rdxd>& Ar2,
NodeValue<Rd>& br1,
NodeValue<Rd>& br2,
const std::vector<NodeId>& map1,
const std::vector<NodeId>& map2) const
{
size_t n = map1.size();
for(size_t i = 0; i<n; i++){
NodeId node_id1 = map1[i]; NodeId node_id2 = map2[i];
Ar1[node_id1] += Ar2[node_id2];
Ar2[node_id2] = Ar1[node_id1];
br1[node_id1] += br2[node_id2];
br2[node_id2] = br1[node_id1];
}
}
template <size_t Dimension>
void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const MeshType& mesh,
NodeValue<Rd>& ur,
NodeValue<Rd>& CR_ur,
NodeValuePerCell<Rd>& Fjr,
NodeValuePerCell<Rd>& CR_Fjr,
const std::vector<NodeId>& map2) const
{
const size_t& n = map2.size();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
for(size_t i = 0; i<n; i++){
const NodeId node_id = map2[i];
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
for(size_t j = 0; j < node_to_cell.size(); j++){
const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
Fjr(J,R) = CR_Fjr(J,R);
}
ur[node_id] = CR_ur[node_id];
}
}
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::FixedBoundaryCondition
{
private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
: m_mesh_node_boundary{mesh_node_boundary}
{}
~FixedBoundaryCondition() = default;
};
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::PressureBoundaryCondition
{
private:
const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
const Array<const double> m_value_list;
public:
const Array<const FaceId>&
faceList() const
{ return m_mesh_face_boundary.faceList();
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
const Array<const double>& value_list)
: m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
{}
~PressureBoundaryCondition() = default;
};
template <>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<1>::PressureBoundaryCondition
{
private:
const MeshNodeBoundary<1> m_mesh_node_boundary;
const Array<const double> m_value_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
: m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
{}
~PressureBoundaryCondition() = default;
};
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition
{
private:
const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
const Array<TinyVector<Dimension>> m_value_list;
const std::shared_ptr<const Socket> m_socket;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
Array<const TinyVector<Dimension>>
valueList(const DiscreteScalarFunction& p) const
{
if (parallel::size() > 1) {
throw NotImplementedError("Parallelism is not supported");
}
if (m_value_list.size() != 1) {
throw NotImplementedError("Non connected boundaries are not supported");
}
const CellId cell_id = m_node_to_cell_matrix[m_mesh_node_boundary.nodeList()[0]][0];
write(*m_socket, p[cell_id]);
read(*m_socket, m_value_list[0]);
return m_value_list;
}
ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh,
const MeshNodeBoundary<Dimension>& mesh_node_boundary,
const std::shared_ptr<const Socket>& socket)
: m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
m_mesh_node_boundary{mesh_node_boundary},
m_value_list(mesh_node_boundary.nodeList().size()),
m_socket{socket}
{
if constexpr (Dimension > 1) {
throw NotImplementedError("external FSI velocity bc in dimension > 1");
}
}
~ExternalFSIVelocityBoundaryCondition() = default;
};
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::VelocityBoundaryCondition
{
private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
const Array<const TinyVector<Dimension>> m_value_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const TinyVector<Dimension>>&
valueList() const
{
return m_value_list;
}
VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
const Array<const TinyVector<Dimension>>& value_list)
: m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
{}
~VelocityBoundaryCondition() = default;
};
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const {
return m_mesh_flat_node_boundary.nodeList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::CouplingBoundaryCondition
{
private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
const size_t m_label;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
size_t
label() const {
return m_label;
}
CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,const size_t label)
: m_mesh_node_boundary{mesh_node_boundary}, m_label{label}
{
;
}
~CouplingBoundaryCondition() = default;
};
LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh1, const std::shared_ptr<const IMesh>& i_mesh2)
{
if (not i_mesh1) {
throw NormalError("mesh1 not defined");
}
if (not i_mesh2) {
throw NormalError("mesh2 not defined");
}
if (not (i_mesh1->dimension()==i_mesh2->dimension())) {
throw NormalError("mesh1 and mesh2 must have the same dimension");
}
switch (i_mesh1->dimension()) {
case 1: {
m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<1>>();
break;
}
case 2: {
m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<2>>(); break;
}
case 3: {
m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<3>>();
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}