diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 049584c707909560d0f9b9e539b8c7c55b919e8a..ece058652cab9c7b95924506e5534f02731bc08b 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -22,6 +22,7 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters)
{
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
@@ -37,21 +38,21 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
{
return createDirichletBoundaryCondition(
config, dof_table, mesh, variable_id,
- integration_order, parameters,
+ integration_order, shapefunction_order, parameters,
mesh_node_searcher, boundary_element_searcher);
}
else if (type == "Neumann")
{
return createNeumannBoundaryCondition(
config, dof_table, mesh, variable_id,
- integration_order, parameters,
+ integration_order, shapefunction_order, parameters,
mesh_node_searcher, boundary_element_searcher);
}
else if (type == "Robin")
{
return createRobinBoundaryCondition(
config, dof_table, mesh, variable_id,
- integration_order, parameters,
+ integration_order, shapefunction_order, parameters,
mesh_node_searcher, boundary_element_searcher);
}
else
@@ -65,6 +66,7 @@ BoundaryConditionBuilder::createDirichletBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned /*integration_order*/,
+ const unsigned /*shapefunction_order*/,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
MeshGeoToolsLib::BoundaryElementsSearcher& /*boundary_element_searcher*/)
@@ -110,6 +112,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
MeshGeoToolsLib::MeshNodeSearcher& /*mesh_node_searcher*/,
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher)
@@ -118,7 +121,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
config.config,
getClonedElements(boundary_element_searcher, config.geometry),
dof_table, variable_id, config.component_id,
- mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+ mesh.isAxiallySymmetric(), integration_order, shapefunction_order, mesh.getDimension(),
parameters);
}
@@ -127,6 +130,7 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
MeshGeoToolsLib::MeshNodeSearcher& /*mesh_node_searcher*/,
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher)
@@ -135,7 +139,7 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
config.config,
getClonedElements(boundary_element_searcher, config.geometry),
dof_table, variable_id, config.component_id,
- mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+ mesh.isAxiallySymmetric(), integration_order, shapefunction_order, mesh.getDimension(),
parameters);
}
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 32f8df74df5e801eb6e9a277a6cd15e95651cdb8..3a6b78039059f82fb8f0a0404423344501bf0cda 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -76,6 +76,7 @@ public:
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters);
protected:
@@ -84,6 +85,7 @@ protected:
const NumLib::LocalToGlobalIndexMap& dof_table,
const MeshLib::Mesh& mesh, const int variable_id,
const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters,
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
@@ -94,6 +96,7 @@ protected:
const NumLib::LocalToGlobalIndexMap& dof_table,
const MeshLib::Mesh& mesh, const int variable_id,
const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters,
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
@@ -104,6 +107,7 @@ protected:
const NumLib::LocalToGlobalIndexMap& dof_table,
const MeshLib::Mesh& mesh, const int variable_id,
const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters,
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
index 4a776f8db17736f509cb0628b2286355402d04d0..b233ed3bc6dd90bee97e1d2213d0abb360198a5d 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
@@ -45,6 +45,8 @@ void DirichletBoundaryCondition::getEssentialBCValues(
// TODO: that might be slow, but only done once
const auto g_idx =
_dof_table.getGlobalIndex(l, _variable_id, _component_id);
+ if (g_idx == NumLib::MeshComponentMap::nop)
+ continue;
// For the DDC approach (e.g. with PETSc option), the negative
// index of g_idx means that the entry by that index is a ghost one,
// which should be dropped. Especially for PETSc routines MatZeroRows
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 515b653b8592c1332358f0e91c71e93cc1e2a89b..25e8988090a56fc94d5194228ff59f66ab5508a3 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -26,6 +26,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
typename std::decay<Data>::type>::value,
bool>::type is_axially_symmetric,
unsigned const integration_order,
+ unsigned const shapefunction_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
@@ -55,7 +56,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
variable_id, component_id, std::move(all_mesh_subsets), _elements));
createLocalAssemblers<LocalAssemblerImplementation>(
- global_dim, _elements, *_dof_table_boundary, _local_assemblers,
+ global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
is_axially_symmetric, _integration_order, _data);
}
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 5edcf191a7302fc8a42439aa265a22fa8700a902..601ec399fcd15c70e8c58ce66d2ac7a141cac3c5 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -33,6 +33,7 @@ public:
typename std::decay<Data>::type>::value,
bool>::type is_axially_symmetric,
unsigned const integration_order,
+ unsigned const shapefunction_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
index 79601405d89f4dc9446a89cdb636d37a45a83a74..7a4f13dfe295afcb166327c7c9e1d12afe519d00 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -17,7 +17,9 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{
DBUG("Constructing Neumann BC from config.");
@@ -31,7 +33,7 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
auto const& param = findParameter<double>(param_name, parameters, 1);
return std::unique_ptr<NeumannBoundaryCondition>(
- new NeumannBoundaryCondition(is_axially_symmetric, integration_order,
+ new NeumannBoundaryCondition(is_axially_symmetric, integration_order, shapefunction_order,
dof_table, variable_id, component_id,
global_dim, std::move(elements), param));
}
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
index 0795935b003ab38a154fdc424bd976b9b8bf7019..c956d824a4ae6f12387895172e2b28ec39c54b87 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
@@ -24,7 +24,9 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters);
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
index eab6bc696704388227a495355aa1148046a88ce0..6a3b9e40d168433dc165b1200b885e41a4f6ccbb 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
@@ -17,7 +17,9 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{
DBUG("Constructing RobinBcConfig from config.");
@@ -33,7 +35,7 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
auto const& u_0 = findParameter<double>(u_0_name, parameters, 1);
return std::unique_ptr<RobinBoundaryCondition>(new RobinBoundaryCondition(
- is_axially_symmetric, integration_order, dof_table, variable_id,
+ is_axially_symmetric, integration_order, shapefunction_order, dof_table, variable_id,
component_id, global_dim, std::move(elements),
RobinBoundaryConditionData{alpha, u_0}));
}
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
index 5c9357f5df290709605dae339e0be8f90d144e3c..64d8f760bd0d4cfe0561aaee31f03abb525ad1c6 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
@@ -40,7 +40,9 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters);
} // ProcessLib
diff --git a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
index 648a448856b7f1d41b0061f346086691b44d4f90..059cd3373cc35ffe6318aad19faaab6fc65f59f7 100644
--- a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
+++ b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
@@ -57,7 +57,7 @@ CalculateSurfaceFlux::CalculateSurfaceFlux(
ProcessLib::createLocalAssemblers<CalculateSurfaceFluxLocalAssembler>(
boundary_mesh.getDimension() + 1, // or bulk_mesh.getDimension()?
- boundary_mesh.getElements(), *dof_table, _local_assemblers,
+ boundary_mesh.getElements(), *dof_table, 1, _local_assemblers,
boundary_mesh.isAxiallySymmetric(), integration_order,
*bulk_element_ids, *bulk_face_ids);
}
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 4170854f7dcbb24703afb80a3fe44fa26777f80b..4bb40f2c8781b0b6442a2435196184d73abd0224 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -44,8 +44,10 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index d481370cf128065b739eb0834b38153e879c58ad..617b90ff0ca5d6804205c92587cc09e380e1bfbb 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -38,8 +38,10 @@ void HeatConductionProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 603d30d776f173bae0aeacfd59fc543cb988d7a6..e06ff0db8b1c4ac87d1de67aaf44c293d57b9d39 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -55,8 +55,10 @@ void LiquidFlowProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
ProcessLib::createLocalAssemblers<LiquidFlowLocalAssembler>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
_gravitational_acceleration, _material_properties);
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index a53c7d6d2f946aefbe7138630b9ba0296fe660df..148e9859d4a68302149285b75e67c25aa35edd4e 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -27,6 +27,8 @@ ProcessVariable::ProcessVariable(
_mesh(mesh),
//! \ogs_file_param{prj__process_variables__process_variable__components}
_n_components(config.getConfigParameter<int>("components")),
+ //! \ogs_file_param{prj__process_variables__process_variable__order}
+ _shapefunction_order(config.getConfigParameter<unsigned>("order")),
_initial_condition(findParameter<double>(
//! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
config.getConfigParameter<std::string>("initial_condition"),
@@ -35,6 +37,9 @@ ProcessVariable::ProcessVariable(
{
DBUG("Constructing process variable %s", _name.c_str());
+ if (_shapefunction_order < 1 || 2 < _shapefunction_order)
+ OGS_FATAL("The given shape function order %d is not supported", _shapefunction_order);
+
// Boundary conditions
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions}
if (auto bcs_config = config.getConfigSubtreeOptional("boundary_conditions"))
@@ -95,6 +100,7 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other)
: _name(std::move(other._name)),
_mesh(other._mesh),
_n_components(other._n_components),
+ _shapefunction_order(other._shapefunction_order),
_initial_condition(std::move(other._initial_condition)),
_bc_configs(std::move(other._bc_configs)),
_bc_builder(std::move(other._bc_builder))
@@ -142,7 +148,7 @@ ProcessVariable::createBoundaryConditions(
for (auto& config : _bc_configs)
bcs.emplace_back(_bc_builder->createBoundaryCondition(
- config, dof_table, _mesh, variable_id, integration_order, parameters));
+ config, dof_table, _mesh, variable_id, integration_order, _shapefunction_order, parameters));
return bcs;
}
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index db444a2fa60d43b51f1a998224fad1c7e1156a1d..b4c2e6570dc22046a76e6e6f7dc898f1d24e1de1 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -64,10 +64,14 @@ public:
// components.
MeshLib::PropertyVector<double>& getOrCreateMeshProperty();
+ unsigned getShapeFunctionOrder() const { return _shapefunction_order; }
+
private:
std::string const _name;
MeshLib::Mesh& _mesh;
const int _n_components;
+ unsigned _shapefunction_order; ///< Order of the shapefunctions. Requires
+ /// appropriate mesh.
Parameter<double> const& _initial_condition;
std::vector<BoundaryConditionConfig> _bc_configs;
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 4a3ce42a33e2d37df624ead517070a4cd8f118f0..6d3b36fd7539a65ab458b9fbb17e012904dd3a4a 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -38,8 +38,10 @@ void RichardsFlowProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
index 9ed876922d76c25872f7f918500c77bbb6893c4f..5cfa0fb5d79bd257bcab583a98d0052435d4e29d 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
@@ -23,6 +23,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
const BoundaryConditionConfig& config,
const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
MeshGeoToolsLib::MeshNodeSearcher& /*mesh_node_searcher*/,
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher)
@@ -31,7 +32,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
config.config,
getClonedElements(boundary_element_searcher, config.geometry),
dof_table, variable_id, config.component_id,
- mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+ mesh.isAxiallySymmetric(), integration_order, shapefunction_order, mesh.getDimension(),
parameters, _fracture_prop);
}
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
index 20fe4f15fde80725247604c25c87f013c653c91f..5e7bac36b86e2bce00a1ef14850703f6f6c24286 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
@@ -45,6 +45,7 @@ private:
const NumLib::LocalToGlobalIndexMap& dof_table,
const MeshLib::Mesh& mesh, const int variable_id,
const unsigned integration_order,
+ const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters,
MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 2eaa6b10edb31225a8924c33ff41b9dad683244b..84222f6cd0b8a0b125fd249fe5eaa0b8627f122d 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -33,6 +33,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
typename std::decay<Data>::type>::value,
bool>::type is_axially_symmetric,
unsigned const integration_order,
+ unsigned const shapefunction_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
unsigned const global_dim,
@@ -63,7 +64,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
variable_id, component_id, std::move(all_mesh_subsets), _elements));
createLocalAssemblers<LocalAssemblerImplementation>(
- global_dim, _elements, *_dof_table_boundary, _local_assemblers,
+ global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
is_axially_symmetric, _integration_order, _data, fracture_prop, variable_id);
}
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 3144eb2b31520aa4cc24cf8e69b8f8b53b772a1b..4a2e677eb794b698a5cbd9fd82c7a0ba3fb5e63b 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -37,6 +37,7 @@ public:
typename std::decay<Data>::type>::value,
bool>::type is_axially_symmetric,
unsigned const integration_order,
+ unsigned const shapefunction_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
unsigned const global_dim,
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
index 5f99aedf0147c6b104958a94b005f8eb3767e50f..fb09e7c3538fd4e9bd1c6aa99482ccebc3ef57f6 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -22,13 +22,13 @@ namespace SmallDeformationWithLIE
using NeumannBoundaryCondition = GenericNaturalBoundaryCondition<
Parameter<double> const&, NeumannBoundaryConditionLocalAssembler>;
-std::unique_ptr<BoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<BoundaryCondition> createNeumannBoundaryCondition(
BaseLib::ConfigTree const& config,
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
FractureProperty const& fracture_prop)
{
@@ -44,7 +44,7 @@ createNeumannBoundaryCondition(
return std::unique_ptr<BoundaryCondition>(
new NeumannBoundaryCondition(
- is_axially_symmetric, integration_order,
+ is_axially_symmetric, integration_order, shapefunction_order,
dof_table, variable_id, component_id,
global_dim, std::move(elements), param, fracture_prop));
}
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
index 5c846581e8b89df312c288e9e2bbf1a8a1f74d41..dc76409af349e8d4a905d7c52f0da3a282eff17f 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
@@ -32,7 +32,9 @@ createNeumannBoundaryCondition(
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
int const component_id, bool is_axially_symmetric,
- unsigned const integration_order, unsigned const global_dim,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
FractureProperty const& fracture_prop);
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 0b39f5cfe275b11e2ed91520a9b84e98e3be97d9..3cef80382c159f6a568b979bed48671ac626c18c 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -140,8 +140,10 @@ void TESProcess::initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh, unsigned const integration_order)
{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
ProcessLib::createLocalAssemblers<TESLocalAssembler>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _assembly_params);
initializeSecondaryVariables();
diff --git a/ProcessLib/Utils/CreateLocalAssemblers.h b/ProcessLib/Utils/CreateLocalAssemblers.h
index 457137ffbebe8c9b0ddb7ec67df4b2e3fb0bcdd8..1dd58de61bd32b123e907cdbddf6ffe60db1f6a2 100644
--- a/ProcessLib/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/Utils/CreateLocalAssemblers.h
@@ -28,6 +28,7 @@ template <unsigned GlobalDim, template <typename, typename, unsigned>
typename LocalAssemblerInterface, typename... ExtraCtorArgs>
void createLocalAssemblers(
NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order,
std::vector<MeshLib::Element*> const& mesh_elements,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args)
@@ -42,7 +43,7 @@ void createLocalAssemblers(
// Populate the vector of local assemblers.
local_assemblers.resize(mesh_elements.size());
- LocalDataInitializer initializer(dof_table);
+ LocalDataInitializer initializer(dof_table, shapefunction_order);
DBUG("Calling local assembler builder for all mesh elements.");
GlobalExecutor::transformDereferenced(
@@ -70,6 +71,7 @@ void createLocalAssemblers(
const unsigned dimension,
std::vector<MeshLib::Element*> const& mesh_elements,
NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args)
{
@@ -79,17 +81,17 @@ void createLocalAssemblers(
{
case 1:
detail::createLocalAssemblers<1, LocalAssemblerImplementation>(
- dof_table, mesh_elements, local_assemblers,
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
case 2:
detail::createLocalAssemblers<2, LocalAssemblerImplementation>(
- dof_table, mesh_elements, local_assemblers,
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
case 3:
detail::createLocalAssemblers<3, LocalAssemblerImplementation>(
- dof_table, mesh_elements, local_assemblers,
+ dof_table, shapefunction_order, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
default:
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
index 944a5f274235c07a75bd71a3f1b540681718bcf1..2708ecc8b479182a2ff45bb6eac6e4adafdcc374 100644
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ b/ProcessLib/Utils/LocalDataInitializer.h
@@ -128,10 +128,16 @@ class LocalDataInitializer final
public:
using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
- explicit LocalDataInitializer(
- NumLib::LocalToGlobalIndexMap const& dof_table)
+ LocalDataInitializer(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ const unsigned shapefunction_order)
: _dof_table(dof_table)
{
+ if (shapefunction_order < 1 || 2 < shapefunction_order)
+ OGS_FATAL("The given shape function order %d is not supported", shapefunction_order);
+
+ if (shapefunction_order == 1)
+ {
// /// Lines and points ///////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
@@ -149,10 +155,9 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Line3))] =
- makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
+ makeLocalAssemblerBuilder<NumLib::ShapeLine2>();
#endif
-
// /// Quads and Hexahedra ///////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \
@@ -170,18 +175,17 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Quad8))] =
- makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
_builder[std::type_index(typeid(MeshLib::Quad9))] =
- makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Hex20))] =
- makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+ makeLocalAssemblerBuilder<NumLib::ShapeHex8>();
#endif
-
// /// Simplices ////////////////////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \
@@ -199,16 +203,15 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Tri6))] =
- makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+ makeLocalAssemblerBuilder<NumLib::ShapeTri3>();
#endif
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Tet10))] =
- makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+ makeLocalAssemblerBuilder<NumLib::ShapeTet4>();
#endif
-
// /// Prisms ////////////////////////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \
@@ -220,7 +223,7 @@ public:
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Prism15))] =
- makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+ makeLocalAssemblerBuilder<NumLib::ShapePrism6>();
#endif
// /// Pyramids //////////////////////////////////////////////////
@@ -231,11 +234,72 @@ public:
makeLocalAssemblerBuilder<NumLib::ShapePyra5>();
#endif
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra5>();
+#endif
+
+ }
+ else if (shapefunction_order == 2)
+ {
+ // /// Lines and points ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_LINE) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 1 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Line3))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
+#endif
+
+
+ // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Quad8))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+ _builder[std::type_index(typeid(MeshLib::Quad9))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Hex20))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+
+ // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tri6))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tet10))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+
+ // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \
+ && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Prism15))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+ // /// Pyramids //////////////////////////////////////////////////
+
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \
&& OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Pyramid13))] =
makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
#endif
+ }
}
/// Sets the provided \c data_ptr to the newly created local assembler data.
diff --git a/Tests/Data b/Tests/Data
index 902f208d67938f094e67b420d0ff994ef8e33957..95d086de5e3c6db78c07002cec395809a157ad3b 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 902f208d67938f094e67b420d0ff994ef8e33957
+Subproject commit 95d086de5e3c6db78c07002cec395809a157ad3b
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index 0189abb0482e09be6d55c50bc964e465dea4bfbc..8ff442621fbd9f22a691f8c9da38014727bf659e 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -146,7 +146,7 @@ public:
// createAssemblers(mesh);
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), *_dof_table,
+ mesh.getDimension(), mesh.getElements(), *_dof_table, 1,
_local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
}