diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 7df72ade933151b4d8c889cb1945b619154c82a8..f92ee5e6f83a5c2fe43f4fad4f692f7107a58815 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -60,18 +60,29 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
project_directory);
detail::readGeometry(geometry_file, *_geoObjects);
- std::string const mesh_file = BaseLib::copyPathToFileName(
+ {
//! \ogs_file_param{prj__mesh}
- project_config.getConfigParameter<std::string>("mesh"),
- project_directory);
+ auto const mesh_param = project_config.getConfigParameter("mesh");
- MeshLib::Mesh* const mesh = MeshLib::IO::readMeshFromFile(mesh_file);
- if (!mesh)
- {
- OGS_FATAL("Could not read mesh from \'%s\' file. No mesh added.",
- mesh_file.c_str());
+ std::string const mesh_file = BaseLib::copyPathToFileName(
+ //! \ogs_file_param{prj__mesh}
+ mesh_param.getValue<std::string>(), project_directory);
+
+ MeshLib::Mesh* const mesh = MeshLib::IO::readMeshFromFile(mesh_file);
+ if (!mesh)
+ {
+ OGS_FATAL("Could not read mesh from \'%s\' file. No mesh added.",
+ mesh_file.c_str());
+ }
+
+ if (auto const axially_symmetric =
+ //! \ogs_file_param{prj__mesh__axial_symmetric}
+ mesh_param.getConfigAttributeOptional<bool>("axially_symmetric"))
+ {
+ mesh->setAxiallySymmetric(*axially_symmetric);
+ }
+ _mesh_vec.push_back(mesh);
}
- _mesh_vec.push_back(mesh);
//! \ogs_file_param{prj__curves}
parseCurves(project_config.getConfigSubtreeOptional("curves"));
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 086222ed57d2e50366a81e84cb459280b2864638..6e275d817f91211a550b298f180e90466f3197af 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -416,6 +416,110 @@ if(NOT OGS_USE_MPI)
cube_1e3_expected_pcs_0_ts_101_t_1.000000.vtu cube_1e3_pcs_0_ts_101_t_1.000000.vtu displacement displacement
)
+ # SQUARE 1x1 GROUNDWATER FLOW TEST -- AXIALLY SYMMETRIC
+ # test results are compared to 3D simulation on a wedge-shaped domain
+ AddTest(
+ NAME GroundWaterFlowProcess_square_1x1_1e2_axi
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_1e2_axi.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1.6e-5 RELTOL 1e-5
+ DIFF_DATA
+ wedge-1e2-ang-0.02-surface.vtu square_1e2_axi_pcs_0_ts_1_t_1.000000.vtu temperature temperature
+ )
+ # # WEDGE 1x1 GROUNDWATER FLOW TEST -- computes reference results for the above test
+ # AddTest(
+ # NAME GroundWaterFlowProcess_wedge_1e2_ang_0.02
+ # PATH Elliptic/square_1x1_GroundWaterFlow
+ # EXECUTABLE ogs
+ # EXECUTABLE_ARGS wedge_1e2_axi_ang_0.02.prj
+ # )
+
+ # SQUARE 1x1 GROUNDWATER FLOW TEST -- AXIALLY SYMMETRIC
+ # test results are compared to 3D simulation on a wedge-shaped domain
+ AddTest(
+ NAME GroundWaterFlowProcess_square_1x1_1e4_axi_ang_0.02
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_1e4_axi.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1.6e-5 RELTOL 1e-5
+ DIFF_DATA
+ wedge-1e4-ang-0.02-surface.vtu square_1e4_axi_pcs_0_ts_1_t_1.000000.vtu temperature temperature
+ )
+ # # WEDGE 1x1 GROUNDWATER FLOW TEST -- computes reference results for the above test
+ # AddTest(
+ # NAME GroundWaterFlowProcess_wedge_1e4_ang_0.02
+ # PATH Elliptic/square_1x1_GroundWaterFlow
+ # EXECUTABLE ogs
+ # EXECUTABLE_ARGS wedge_1e4_axi_ang_0.02.prj
+ # )
+
+ # SMALL DEFORMATION TEST -- AXIALLY SYMMETRIC
+ AddTest(
+ NAME SmallDeformation_ring_plane_strain_axi
+ PATH Mechanics/Linear
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS ring_plane_strain.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 6e-4 RELTOL 1e-4
+ DIFF_DATA
+ ring_plane_strain_1e4_solution.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu u displacement
+ ring_plane_strain_1e4_solution.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu sigma_rr sigma_rr
+ ring_plane_strain_1e4_solution.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu sigma_ff sigma_ff
+ )
+
+ # SQUARE 1x1 HEAT CONDUCTION TEST -- AXIALLY SYMMETRIC
+ # test results are compared to 3D simulation on a wedge-shaped domain
+ AddTest(
+ NAME 2D_HeatConduction_axi
+ PATH Parabolic/T/2D_axially_symmetric
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_1e2_axi.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1.7e-5 RELTOL 1e-5
+ DIFF_DATA
+ wedge_1e2_axi_ang_0.02_t_2s_extracted_surface.vtu square_1e2_axi_pcs_0_ts_2_t_2.000000.vtu temperature temperature
+ wedge_1e2_axi_ang_0.02_t_2s_extracted_surface.vtu square_1e2_axi_pcs_0_ts_2_t_2.000000.vtu heat_flux_x heat_flux_x
+ )
+ # # WEDGE 1x1 HEATCONDUCTION TEST -- computes reference results for the above test
+ # AddTest(
+ # NAME 2D_HeatConduction_wedge
+ # PATH Parabolic/T/2D_axially_symmetric
+ # EXECUTABLE ogs
+ # EXECUTABLE_ARGS wedge_1e2_axi_ang_0.02.prj
+ # )
+
+ # SQUARE 1x1 TES TEST -- AXIALLY SYMMETRIC
+ # test results are compared to 3D simulation on a wedge-shaped domain
+ AddTest(
+ NAME LARGE_TES_inert_axi
+ PATH Parabolic/TES/2D
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS tes-inert-axi.prj
+ WRAPPER time
+ TESTER vtkdiff
+ # Note: Since the temperature and pressure only vary by a factor of ~ 1.e-6 in x-direction
+ # the relative tolerance has to be much smaller than 1.e-6
+ ABSTOL 1e-12 RELTOL 2e-9
+ DIFF_DATA
+ inert-wedge-extracted-surface-t-1s.vtu tes_inert_axi_pcs_0_ts_4_t_1.000000.vtu pressure pressure
+ inert-wedge-extracted-surface-t-1s.vtu tes_inert_axi_pcs_0_ts_4_t_1.000000.vtu temperature temperature
+ inert-wedge-extracted-surface-t-1s.vtu tes_inert_axi_pcs_0_ts_4_t_1.000000.vtu v_mass_frac v_mass_frac
+ )
+ # # WEDGE 1x1 TES TEST -- computes reference results for the above test
+ # AddTest(
+ # NAME TES_inert_wedge
+ # PATH Parabolic/TES/2D
+ # EXECUTABLE ogs
+ # EXECUTABLE_ARGS tes-inert-wedge.prj
+ # )
+
else()
# MPI groundwater flow tests
AddTest(
diff --git a/CHANGELOG.md b/CHANGELOG.md
index eccdfb5b4107c557900b3306570aff964d666b1c..91f206929369d7bb7ca957859d8d9666bc16e976 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -26,6 +26,8 @@ for boundary conditions were generalized.
ID dependent values. #1426
- Added fluid property class and several fluid density and viscosity models
based on it. #1398, #1435
+ - Enabled solving of axially symmetric problems on 2D meshes for all currently
+ implemented processes. #1443
#### Utilities
New utilities:
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index 4b7a72bf0adc27b05ded75b4c677337ba73f058b..0a4e215901043e5b46269195d24e70a0dc33c37c 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -125,6 +125,16 @@ public:
MeshLib::Properties & getProperties() { return _properties; }
MeshLib::Properties const& getProperties() const { return _properties; }
+ bool isAxiallySymmetric() const { return _is_axially_symmetric; }
+ void setAxiallySymmetric(bool is_axial_symmetric) {
+ _is_axially_symmetric = is_axial_symmetric;
+ if (_is_axially_symmetric && getDimension() != 2) {
+ OGS_FATAL(
+ "Axial symmetry is implemented only for two-dimensional "
+ "meshes.");
+ }
+ }
+
protected:
/// Set the minimum and maximum length over the edges of the mesh.
void calcEdgeLengthRange();
@@ -163,6 +173,8 @@ protected:
std::vector<Element*> _elements;
std::size_t _n_base_nodes;
Properties _properties;
+
+ bool _is_axially_symmetric = false;
}; /* class */
diff --git a/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h b/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
index eeb5ee23787a6e716fe830ace6048e1bd2186bf1..89f010bf63e536dd163bb8b6d88abd98e34c2ac4 100644
--- a/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
+++ b/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
@@ -64,6 +64,7 @@ inline void setZero(ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX> &shape, ShapeDataFie
setZero(shape, ShapeDataFieldType<ShapeMatrixType::DNDR>());
setMatrixZero(shape.J);
shape.detJ = .0;
+ shape.integralMeasure = 0.0;
}
template <class T_N, class T_DNDR, class T_J, class T_DNDX>
diff --git a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
index a1d6d5d8ca46982ce12ea9d1e58a711f0ef3037e..c1c60d3730f0c8a8f0b390c687355a92d05beb2d 100644
--- a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
+++ b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
@@ -54,6 +54,7 @@ struct ShapeMatrices
double detJ; ///< Determinant of the Jacobian
JacobianType invJ; ///< Inverse matrix of the Jacobian
DxShapeType dNdx; ///< Matrix of gradient of shape functions in physical coordinates, dN(r)/dx
+ double integralMeasure;
/** Not default constructible, dimensions always must be given.
*
@@ -80,7 +81,8 @@ struct ShapeMatrices
J(local_dim, local_dim),
detJ(.0),
invJ(local_dim, local_dim),
- dNdx(global_dim, n_nodes)
+ dNdx(global_dim, n_nodes),
+ integralMeasure(0.0)
{
setZero();
}
diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index df6de49cb2875e8765645f5cf622482f8ba346cd..6f804a7bb74c378b46e6251e038e5b6ed6bece43 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -16,9 +16,10 @@
#include <cassert>
+#include <boost/math/constants/constants.hpp>
-#include "../CoordinatesMapping/ShapeMatrices.h"
-#include "../CoordinatesMapping/NaturalCoordinatesMapping.h"
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
+#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
namespace NumLib
{
@@ -85,10 +86,12 @@ public:
* @param global_dim global dimension
*/
void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
- const unsigned global_dim) const
+ const unsigned global_dim,
+ bool is_axially_symmetric) const
{
NaturalCoordsMappingType::computeShapeMatrices(*_ele, natural_pt, shape,
global_dim);
+ computeIntegralMeasure(is_axially_symmetric, shape);
}
/**
@@ -101,13 +104,46 @@ public:
*/
template <ShapeMatrixType T_SHAPE_MATRIX_TYPE>
void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
- const unsigned global_dim) const
+ const unsigned global_dim,
+ bool is_axially_symmetric) const
{
NaturalCoordsMappingType::template computeShapeMatrices<
T_SHAPE_MATRIX_TYPE>(*_ele, natural_pt, shape, global_dim);
+ computeIntegralMeasure(is_axially_symmetric, shape);
+ }
+
+ double interpolateZerothCoordinate(
+ typename ShapeMatrices::ShapeType const& N) const
+ {
+ auto* nodes = _ele->getNodes();
+ typename ShapeMatrices::ShapeType rs(N.size());
+ for (int i=0; i<rs.size(); ++i) {
+ rs[i] = (*nodes[i])[0];
+ }
+ auto const r = N.dot(rs);
+ return r;
}
private:
+ void computeIntegralMeasure(bool is_axially_symmetric,
+ ShapeMatrices& shape) const
+ {
+ if (!is_axially_symmetric) {
+ shape.integralMeasure = 1.0;
+ return;
+ }
+
+ // Note: If an integration point is located at the rotation axis, r will
+ // be zero which might lead to problems with the assembled equation
+ // system.
+ // E.g., for triangle elements, if an integration point is
+ // located at edge of the unit triangle, it is possible that
+ // r becomes zero.
+ auto const r = interpolateZerothCoordinate(shape.N);
+ shape.integralMeasure =
+ boost::math::constants::two_pi<double>() * r;
+ }
+
const MeshElementType* _ele;
};
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 311377ea6df090f0d537575ac2b1eba796fea326..fbe44992065514d81e342c683ce80d6d9c1ba088 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -85,15 +85,17 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
return createNeumannBoundaryCondition(
config.config,
getClonedElements(boundary_element_searcher, config.geometry),
- dof_table, variable_id, config.component_id, integration_order,
- mesh.getDimension(), parameters);
+ dof_table, variable_id, config.component_id,
+ mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+ parameters);
}
else if (type == "Robin") {
return createRobinBoundaryCondition(
config.config,
getClonedElements(boundary_element_searcher, config.geometry),
- dof_table, variable_id, config.component_id, integration_order,
- mesh.getDimension(), parameters);
+ dof_table, variable_id, config.component_id,
+ mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+ parameters);
}
else
{
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 8c0db0201b83c563b2eca8cd537889a8532d09e0..515b653b8592c1332358f0e91c71e93cc1e2a89b 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -24,11 +24,12 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
typename std::enable_if<
std::is_same<typename std::decay<BoundaryConditionData>::type,
typename std::decay<Data>::type>::value,
- unsigned const>::type integration_order,
+ bool>::type is_axially_symmetric,
+ unsigned const integration_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
- unsigned const global_dim,
- std::vector<MeshLib::Element*>&& elements, Data&& data)
+ unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+ Data&& data)
: _data(std::forward<Data>(data)),
_elements(std::move(elements)),
_integration_order(integration_order)
@@ -54,8 +55,8 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
variable_id, component_id, std::move(all_mesh_subsets), _elements));
createLocalAssemblers<LocalAssemblerImplementation>(
- global_dim, _elements, *_dof_table_boundary, _integration_order,
- _local_assemblers, _data);
+ global_dim, _elements, *_dof_table_boundary, _local_assemblers,
+ is_axially_symmetric, _integration_order, _data);
}
template <typename BoundaryConditionData,
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index f40fc809009a51d35509888a8ce32c9948623017..5edcf191a7302fc8a42439aa265a22fa8700a902 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -31,11 +31,12 @@ public:
typename std::enable_if<
std::is_same<typename std::decay<BoundaryConditionData>::type,
typename std::decay<Data>::type>::value,
- unsigned const>::type integration_order,
+ bool>::type is_axially_symmetric,
+ unsigned const integration_order,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id,
- unsigned const global_dim,
- std::vector<MeshLib::Element*>&& elements, Data&& data);
+ unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+ Data&& data);
~GenericNaturalBoundaryCondition() override;
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
index bac5f8e35a495dfcfe98319dbfa54f10bac69d96..523e55abed423607a6d7539239b3b8529f87d924 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
@@ -38,11 +38,12 @@ protected:
public:
GenericNaturalBoundaryConditionLocalAssembler(
- MeshLib::Element const& e, unsigned const integration_order)
+ MeshLib::Element const& e, bool is_axially_symmetric,
+ unsigned const integration_order)
: _integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
- e, _integration_method))
+ e, is_axially_symmetric, _integration_method))
{
}
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
index 448ef07b5354b2ae2a07c555f30b54d8a20ef8f5..79601405d89f4dc9446a89cdb636d37a45a83a74 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -12,13 +12,12 @@
namespace ProcessLib
{
-std::unique_ptr<NeumannBoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
BaseLib::ConfigTree const& config,
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, unsigned const integration_order,
- unsigned const global_dim,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{
DBUG("Constructing Neumann BC from config.");
@@ -32,9 +31,9 @@ createNeumannBoundaryCondition(
auto const& param = findParameter<double>(param_name, parameters, 1);
return std::unique_ptr<NeumannBoundaryCondition>(
- new NeumannBoundaryCondition(
- integration_order, dof_table, variable_id, component_id,
- global_dim, std::move(elements), param));
+ new NeumannBoundaryCondition(is_axially_symmetric, integration_order,
+ dof_table, variable_id, component_id,
+ global_dim, std::move(elements), param));
}
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
index 53713fb255bda0488943968e0978deae81afa799..0795935b003ab38a154fdc424bd976b9b8bf7019 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
@@ -19,13 +19,12 @@ namespace ProcessLib
using NeumannBoundaryCondition = GenericNaturalBoundaryCondition<
Parameter<double> const&, NeumannBoundaryConditionLocalAssembler>;
-std::unique_ptr<NeumannBoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
BaseLib::ConfigTree const& config,
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, unsigned const integration_order,
- unsigned const global_dim,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters);
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index fd9a5d43e9beff7eff9f6e0d5623d366a7e894b1..2348484d577b2e1ad54e6a9c7645bb7c8d6c7f19 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -31,9 +31,10 @@ public:
NeumannBoundaryConditionLocalAssembler(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
unsigned const integration_order,
Parameter<double> const& neumann_bc_parameter)
- : Base(e, integration_order),
+ : Base(e, is_axially_symmetric, integration_order),
_neumann_bc_parameter(neumann_bc_parameter),
_local_rhs(local_matrix_size)
{
@@ -58,7 +59,8 @@ public:
auto const& sm = Base::_shape_matrices[ip];
auto const& wp = Base::_integration_method.getWeightedPoint(ip);
_local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] *
- sm.detJ * wp.getWeight();
+ sm.detJ * wp.getWeight() *
+ sm.integralMeasure;
}
auto const indices = NumLib::getIndices(id, dof_table_boundary);
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
index df8c74183ae35cc60463ae80d0db574643ba5ef1..eab6bc696704388227a495355aa1148046a88ce0 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
@@ -16,8 +16,8 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
BaseLib::ConfigTree const& config,
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, unsigned const integration_order,
- unsigned const global_dim,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters)
{
DBUG("Constructing RobinBcConfig from config.");
@@ -32,11 +32,10 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
auto const& alpha = findParameter<double>(alpha_name, parameters, 1);
auto const& u_0 = findParameter<double>(u_0_name, parameters, 1);
- return std::unique_ptr<RobinBoundaryCondition>(
- new RobinBoundaryCondition(
- integration_order, dof_table, variable_id, component_id, global_dim,
- std::move(elements),
- RobinBoundaryConditionData{alpha, u_0}));
+ return std::unique_ptr<RobinBoundaryCondition>(new RobinBoundaryCondition(
+ is_axially_symmetric, integration_order, dof_table, variable_id,
+ component_id, global_dim, std::move(elements),
+ RobinBoundaryConditionData{alpha, u_0}));
}
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
index 2b67f531c018b70a476c38bc8fe26137dea3431c..5c9357f5df290709605dae339e0be8f90d144e3c 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
@@ -39,8 +39,8 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
BaseLib::ConfigTree const& config,
std::vector<MeshLib::Element*>&& elements,
NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, unsigned const integration_order,
- unsigned const global_dim,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const global_dim,
std::vector<std::unique_ptr<ParameterBase>> const& parameters);
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index a6322168e37758d3a058cdf6684af416ff7f6497..90ab60ceab4b9adc9ab3073f319b45b3a08ded7e 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -31,11 +31,12 @@ class RobinBoundaryConditionLocalAssembler final
ShapeFunction, IntegrationMethod, GlobalDim>;
public:
- RobinBoundaryConditionLocalAssembler(
- MeshLib::Element const& e, std::size_t const local_matrix_size,
- unsigned const integration_order,
- RobinBoundaryConditionData const& data)
- : Base(e, integration_order),
+ RobinBoundaryConditionLocalAssembler(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ RobinBoundaryConditionData const& data)
+ : Base(e, is_axially_symmetric, integration_order),
_data(data),
_local_K(local_matrix_size, local_matrix_size),
_local_rhs(local_matrix_size)
@@ -70,9 +71,9 @@ public:
// adding a alpha term to the diagonal of the stiffness matrix
// and a alpha * u_0 term to the rhs vector
_local_K.diagonal().noalias() +=
- sm.N * alpha * sm.detJ * wp.getWeight();
- _local_rhs.noalias() +=
- sm.N * alpha * u_0 * sm.detJ * wp.getWeight();
+ sm.N * alpha * sm.detJ * wp.getWeight() * sm.integralMeasure;
+ _local_rhs.noalias() += sm.N * alpha * u_0 * sm.detJ *
+ wp.getWeight() * sm.integralMeasure;
}
auto const indices = NumLib::getIndices(id, dof_table_boundary);
diff --git a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
index bd730371cbbeb103f4e8e4a8d0b0ac7b38e81626..51557f5d4ed3e3964324e46c5048414f74025c8b 100644
--- a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
+++ b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
@@ -48,11 +48,12 @@ CalculateSurfaceFlux::CalculateSurfaceFlux(MeshLib::Mesh& boundary_mesh,
boost::optional<MeshLib::PropertyVector<std::size_t> const&> bulk_face_ids(
boundary_mesh.getProperties().template getPropertyVector<std::size_t>(
"OriginalFaceIDs"));
- const std::size_t integration_order = 2;
+ const unsigned integration_order = 2;
ProcessLib::createLocalAssemblers<CalculateSurfaceFluxLocalAssembler>(
- boundary_mesh.getDimension()+1, // or bulk_mesh.getDimension()?
- boundary_mesh.getElements(), *dof_table, integration_order,
- _local_assemblers, *bulk_element_ids, *bulk_face_ids);
+ boundary_mesh.getDimension() + 1, // or bulk_mesh.getDimension()?
+ boundary_mesh.getElements(), *dof_table, _local_assemblers,
+ boundary_mesh.isAxiallySymmetric(), integration_order,
+ *bulk_element_ids, *bulk_face_ids);
}
void CalculateSurfaceFlux::integrate(GlobalVector const& x,
diff --git a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFluxLocalAssembler.h b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFluxLocalAssembler.h
index f06b404a5817f0e039d1da9df1e97a8a7a5469a0..5dbe274ef0d7fc412e8649c8c69af40c194e3545 100644
--- a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFluxLocalAssembler.h
+++ b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFluxLocalAssembler.h
@@ -55,6 +55,7 @@ public:
CalculateSurfaceFluxLocalAssembler(
MeshLib::Element const& surface_element,
std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
unsigned const integration_order,
MeshLib::PropertyVector<std::size_t> const& bulk_element_ids,
MeshLib::PropertyVector<std::size_t> const& bulk_face_ids)
@@ -74,15 +75,16 @@ public:
std::vector<typename ShapeMatricesType::ShapeMatrices> shape_matrices;
shape_matrices.reserve(n_integration_points);
- _detJ.reserve(n_integration_points);
+ _detJ_times_integralMeasure.reserve(n_integration_points);
for (std::size_t ip = 0; ip < n_integration_points; ++ip)
{
shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
ShapeFunction::NPOINTS);
fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N_J>(
_integration_method.getWeightedPoint(ip).getCoords(),
- shape_matrices[ip], GlobalDim);
- _detJ.push_back(shape_matrices[ip].detJ);
+ shape_matrices[ip], GlobalDim, is_axially_symmetric);
+ _detJ_times_integralMeasure.push_back(
+ shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure);
}
}
@@ -133,7 +135,8 @@ public:
surface_element_normal.getCoords(), 3)));
balance.getComponent(element_id, component_id) +=
- bulk_grad_times_normal * _detJ[ip] * wp.getWeight();
+ bulk_grad_times_normal * _detJ_times_integralMeasure[ip] *
+ wp.getWeight();
}
}
}
@@ -141,7 +144,7 @@ public:
private:
MeshLib::Element const& _surface_element;
- std::vector<double> _detJ;
+ std::vector<double> _detJ_times_integralMeasure;
IntegrationMethod const _integration_method;
std::size_t _bulk_element_id;
diff --git a/ProcessLib/Deformation/LinearBMatrix.h b/ProcessLib/Deformation/LinearBMatrix.h
index 67a62e072feca45d176456239b55811f78877639..f7ce63a010b10f3e1326907c9ab5f6b721356bc2 100644
--- a/ProcessLib/Deformation/LinearBMatrix.h
+++ b/ProcessLib/Deformation/LinearBMatrix.h
@@ -16,10 +16,31 @@ namespace ProcessLib
{
namespace LinearBMatrix
{
+namespace detail
+{
+template <int NPOINTS, typename DNDX_Type, typename BMatrixType>
+void fillBMatrix2DCartesianPart(DNDX_Type const& dNdx, BMatrixType& b_matrix)
+{
+ for (int i = 0; i < NPOINTS; ++i) {
+ b_matrix(1, NPOINTS + i) = dNdx(1, i);
+ b_matrix(3, i) = dNdx(1, i) / std::sqrt(2);
+ b_matrix(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
+ b_matrix(0, i) = dNdx(0, i);
+ }
+}
+} // detail
+
/// Fills a B-matrix based on given shape function dN/dx values.
-template <int DisplacementDim, int NPOINTS, typename DNDX_Type,
+template <int DisplacementDim,
+ int NPOINTS,
+ typename N_Type,
+ typename DNDX_Type,
typename BMatrixType>
-void computeBMatrix(DNDX_Type const& dNdx, BMatrixType& b_matrix)
+void computeBMatrix(DNDX_Type const& dNdx,
+ BMatrixType& b_matrix,
+ const bool is_axially_symmetric,
+ N_Type const& N,
+ const double radius)
{
static_assert(0 < DisplacementDim && DisplacementDim <= 3,
"LinearBMatrix::computeBMatrix: DisplacementDim must be in "
@@ -38,14 +59,16 @@ void computeBMatrix(DNDX_Type const& dNdx, BMatrixType& b_matrix)
b_matrix(5, i) = dNdx(2, i) / std::sqrt(2);
b_matrix(5, 2 * NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
}
- // no break for fallthrough.
+ detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, b_matrix);
+ break;
case 2:
- for (int i = 0; i < NPOINTS; ++i)
+ detail::fillBMatrix2DCartesianPart<NPOINTS>(dNdx, b_matrix);
+ if (is_axially_symmetric)
{
- b_matrix(1, NPOINTS + i) = dNdx(1, i);
- b_matrix(3, i) = dNdx(1, i) / std::sqrt(2);
- b_matrix(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
- b_matrix(0, i) = dNdx(0, i);
+ for (int i = 0; i < NPOINTS; ++i)
+ {
+ b_matrix(2, i) = N[i] / radius;
+ }
}
break;
default:
diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
index 12a538c53cb3e1b0bae81851c36e25d56d33371d..8a4cb4fd0ab1a72ceae941ba54e768261433d589 100644
--- a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
@@ -82,6 +82,9 @@ std::unique_ptr<Process> createGroundwaterFlowProcess(
"\"%s\"\n\toutput to: \"%s\"",
mesh_name.c_str(), balance_pv_name.c_str(),
balance_out_fname.c_str());
+
+ // Surface mesh and bulk mesh must have equal axial symmetry flags!
+ surface_mesh->setAxiallySymmetric(mesh.isAxiallySymmetric());
}
return std::unique_ptr<Process>{new GroundwaterFlowProcess{
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 48f747953956a37374a3930d80707e64b40454fe..ab55e14a9d060547ca52801b04356331e078481f 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -66,6 +66,7 @@ public:
/// element matrix.
LocalAssemblerData(MeshLib::Element const& element,
std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
unsigned const integration_order,
GroundwaterFlowProcessData const& process_data)
: _element(element),
@@ -73,7 +74,7 @@ public:
_integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
- element, _integration_method))
+ element, is_axially_symmetric, _integration_method))
{
}
@@ -102,8 +103,8 @@ public:
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
- local_K.noalias() +=
- sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight();
+ local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
+ sm.integralMeasure * wp.getWeight();
// Darcy velocity only computed for output.
GlobalDimVectorType const darcy_velocity =
@@ -134,8 +135,10 @@ public:
typename ShapeMatricesType::ShapeMatrices shape_matrices(
ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+ // Note: Axial symmetry is set to false here, because we only need dNdx
+ // here, which is not affected by axial symmetry.
fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
- GlobalDim);
+ GlobalDim, false);
std::vector<double> flux;
flux.resize(3);
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 7fd6d066dff01e211e6d1ec92d6063dd0a011e20..38c6120716ecdf8cec033bf93f1273a6fb357cf9 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -44,8 +44,8 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
unsigned const integration_order)
{
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
- _local_assemblers, _process_data);
+ mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
"darcy_velocity_x", 1,
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index f4cdd26d1861710a63e1b12ac953ad1258623797..96994733b56c15363c6bfb09db786d6a4e73b875 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -62,6 +62,7 @@ public:
/// element matrix.
LocalAssemblerData(MeshLib::Element const& element,
std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
unsigned const integration_order,
HeatConductionProcessData const& process_data)
: _element(element),
@@ -69,11 +70,12 @@ public:
_integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
- element, _integration_method))
+ element, is_axially_symmetric, _integration_method))
{
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+ (void) local_matrix_size;
}
void assemble(double const t, std::vector<double> const& local_x,
@@ -105,10 +107,11 @@ public:
auto const heat_capacity = _process_data.heat_capacity(t, pos)[0];
auto const density = _process_data.density(t, pos)[0];
- local_K.noalias() +=
- sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight();
+ local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
+ wp.getWeight() * sm.integralMeasure;
local_M.noalias() += sm.N.transpose() * density * heat_capacity *
- sm.N * sm.detJ * wp.getWeight();
+ sm.N * sm.detJ * wp.getWeight() *
+ sm.integralMeasure;
// heat flux only computed for output.
GlobalDimVectorType const heat_flux =
-k * sm.dNdx * Eigen::Map<const NodalVectorType>(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 8523f8648701fe0ccf280066e2295343c6dd902e..4b556ed908e7185d5a2461c512ee1747bec6c12f 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -38,8 +38,8 @@ void HeatConductionProcess::initializeConcreteProcess(
unsigned const integration_order)
{
ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
- _local_assemblers, _process_data);
+ mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
"heat_flux_x", 1,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index f1f887c4f497e90c59480ea01fb2deec13fcd5e6..fc0520fd0005a190e9c8a1ec06fbdcdf490f9322 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -57,8 +57,8 @@ public:
/// Computes the flux in the point \c p_local_coords that is given in local
/// coordinates using the values from \c local_x.
virtual std::vector<double> getFlux(
- MathLib::Point3d const& p_local_coords,
- std::vector<double> const& local_x) const
+ MathLib::Point3d const& /*p_local_coords*/,
+ std::vector<double> const& /*local_x*/) const
{
return std::vector<double>();
}
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 8e48390979f624281242a26d57fcf42f12c7f4e9..7139941fb81778eb368534bfbecfc7583004fb33 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -103,9 +103,9 @@ public:
}
// used as call back for CalculateSurfaceFlux process
- virtual std::vector<double> getFlux(std::size_t element_id,
- MathLib::Point3d const& p,
- GlobalVector const&) const
+ virtual std::vector<double> getFlux(std::size_t /*element_id*/,
+ MathLib::Point3d const& /*p*/,
+ GlobalVector const& /*x*/) const
{
return std::vector<double>{};
}
diff --git a/ProcessLib/SmallDeformation/CreateLocalAssemblers.h b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
index aaa61b6933bd5883048bde59b379c90faa747951..f355329b675012abec35ec09a8467a27c1236401 100644
--- a/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
+++ b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
@@ -30,7 +30,6 @@ template <unsigned GlobalDim, int DisplacementDim,
void createLocalAssemblers(
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<MeshLib::Element*> const& mesh_elements,
- unsigned const integration_order,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args)
{
@@ -51,7 +50,6 @@ void createLocalAssemblers(
initializer,
mesh_elements,
local_assemblers,
- integration_order,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
}
@@ -76,7 +74,6 @@ void createLocalAssemblers(
const unsigned dimension,
std::vector<MeshLib::Element*> const& mesh_elements,
NumLib::LocalToGlobalIndexMap const& dof_table,
- unsigned const integration_order,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args)
{
@@ -87,13 +84,13 @@ void createLocalAssemblers(
case 2:
detail::createLocalAssemblers<2, DisplacementDim,
LocalAssemblerImplementation>(
- dof_table, mesh_elements, integration_order, local_assemblers,
+ dof_table, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
case 3:
detail::createLocalAssemblers<3, DisplacementDim,
LocalAssemblerImplementation>(
- dof_table, mesh_elements, integration_order, local_assemblers,
+ dof_table, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
default:
diff --git a/ProcessLib/SmallDeformation/LocalDataInitializer.h b/ProcessLib/SmallDeformation/LocalDataInitializer.h
index 59363f8e5c302ac49868f95208744dd4495639b9..b17a932ac35f6bdd853cd6b8a7ccbfdf7e80b3f6 100644
--- a/ProcessLib/SmallDeformation/LocalDataInitializer.h
+++ b/ProcessLib/SmallDeformation/LocalDataInitializer.h
@@ -219,7 +219,6 @@ public:
std::size_t const id,
MeshLib::Element const& mesh_item,
LADataIntfPtr& data_ptr,
- unsigned const integration_order,
ConstructorArgs&&... args) const
{
auto const type_idx = std::type_index(typeid(mesh_item));
@@ -228,7 +227,7 @@ public:
if (it != _builder.end()) {
auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
data_ptr = it->second(
- mesh_item, num_local_dof, integration_order,
+ mesh_item, num_local_dof,
std::forward<ConstructorArgs>(args)...);
} else {
OGS_FATAL("You are trying to build a local assembler for an unknown mesh element type (%s)."
@@ -241,7 +240,6 @@ private:
using LADataBuilder = std::function<LADataIntfPtr(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
- unsigned const integration_order,
ConstructorArgs&&...
)>;
@@ -260,12 +258,11 @@ private:
{
return [](MeshLib::Element const& e,
std::size_t const local_matrix_size,
- unsigned const integration_order,
ConstructorArgs&&... args)
{
return LADataIntfPtr{
new LAData<ShapeFct>{
- e, local_matrix_size, integration_order,
+ e, local_matrix_size,
std::forward<ConstructorArgs>(args)...
}};
};
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 01cb3ab818c3fc424ab74420f152ad1622342bed..c44e7fa686baeeccfafa133a6b266ca3042ea370 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -51,7 +51,8 @@ struct IntegrationPointData final
_solid_material(other._solid_material),
_material_state_variables(std::move(other._material_state_variables)),
_C(std::move(other._C)),
- _detJ(std::move(other._detJ))
+ _detJ(std::move(other._detJ)),
+ _integralMeasure(other._integralMeasure)
{
}
#endif // _MSC_VER
@@ -67,6 +68,7 @@ struct IntegrationPointData final
typename BMatricesType::KelvinMatrixType _C;
double _detJ;
+ double _integralMeasure;
void pushBackState()
{
@@ -87,9 +89,6 @@ template <typename ShapeMatrixType>
struct SecondaryData
{
std::vector<ShapeMatrixType> N;
- std::vector<double> _sigmaXX;
- std::vector<double> _sigmaYY;
- std::vector<double> _sigmaXY;
};
struct SmallDeformationLocalAssemblerInterface
@@ -97,13 +96,22 @@ struct SmallDeformationLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
virtual std::vector<double> const& getIntPtSigmaXX(
- std::vector<double>& /*cache*/) const = 0;
+ std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYY(
- std::vector<double>& /*cache*/) const = 0;
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaZZ(
+ std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXY(
- std::vector<double>& /*cache*/) const = 0;
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaXZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaYZ(
+ std::vector<double>& cache) const = 0;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -132,6 +140,7 @@ public:
SmallDeformationLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
unsigned const integration_order,
SmallDeformationProcessData<DisplacementDim>& process_data)
: _process_data(process_data),
@@ -142,29 +151,31 @@ public:
_integration_method.getNumberOfPoints();
_ip_data.reserve(n_integration_points);
-
_secondary_data.N.resize(n_integration_points);
- _secondary_data._sigmaXX.resize(n_integration_points);
- _secondary_data._sigmaYY.resize(n_integration_points);
- _secondary_data._sigmaXY.resize(n_integration_points);
auto const shape_matrices =
initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, DisplacementDim>(
- e, _integration_method);
+ e, is_axially_symmetric, _integration_method);
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(*_process_data._material);
auto& ip_data = _ip_data[ip];
- ip_data._detJ = shape_matrices[ip].detJ;
+ auto const& sm = shape_matrices[ip];
+ ip_data._detJ = sm.detJ;
+ ip_data._integralMeasure = sm.integralMeasure;
ip_data._b_matrices.resize(
KelvinVectorDimensions<DisplacementDim>::value,
ShapeFunction::NPOINTS * DisplacementDim);
- LinearBMatrix::computeBMatrix<
- DisplacementDim, ShapeFunction::NPOINTS,
- typename ShapeMatricesType::GlobalDimNodalMatrixType,
- BMatrixType>(shape_matrices[ip].dNdx, ip_data._b_matrices);
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(e,
+ sm.N);
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS>(
+ sm.dNdx, ip_data._b_matrices, is_axially_symmetric, sm.N,
+ x_coord);
ip_data._sigma.resize(
KelvinVectorDimensions<DisplacementDim>::value);
@@ -218,6 +229,7 @@ public:
x_position.setIntegrationPoint(ip);
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const& detJ = _ip_data[ip]._detJ;
+ auto const& integralMeasure = _ip_data[ip]._integralMeasure;
auto const& B = _ip_data[ip]._b_matrices;
auto const& eps_prev = _ip_data[ip]._eps_prev;
@@ -239,15 +251,10 @@ public:
sigma, C, material_state_variables))
OGS_FATAL("Computation of local constitutive relation failed.");
- local_b.noalias() -= B.transpose() * sigma * detJ * wp.getWeight();
+ local_b.noalias() -=
+ B.transpose() * sigma * detJ * wp.getWeight() * integralMeasure;
local_Jac.noalias() +=
- B.transpose() * C * B * detJ * wp.getWeight();
-
-
- // TODO: Reuse _ip_data[ip]._sigma
- _secondary_data._sigmaXX[ip] = sigma[0];
- _secondary_data._sigmaYY[ip] = sigma[1];
- _secondary_data._sigmaXY[ip] = sigma[3];
+ B.transpose() * C * B * detJ * wp.getWeight() * integralMeasure;
}
}
@@ -274,24 +281,57 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
- std::vector<double>& /*cache*/) const override
+ std::vector<double>& cache) const override
{
- return _secondary_data._sigmaXX;
+ return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
- std::vector<double>& /*cache*/) const override
+ std::vector<double>& cache) const override
{
- return _secondary_data._sigmaYY;
+ return getIntPtSigma(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtSigmaZZ(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
- std::vector<double>& /*cache*/) const override
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtSigmaXZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtSigmaYZ(
+ std::vector<double>& cache) const override
{
- return _secondary_data._sigmaXY;
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 5);
}
private:
+ std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
+ std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data) {
+ cache.push_back(ip_data._sigma[component]);
+ }
+
+ return cache;
+ }
+
SmallDeformationProcessData<DisplacementDim>& _process_data;
std::vector<IntegrationPointData<BMatricesType, DisplacementDim>> _ip_data;
@@ -314,11 +354,13 @@ public:
LocalAssemblerData(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
unsigned const integration_order,
SmallDeformationProcessData<DisplacementDim>& process_data)
: SmallDeformationLocalAssembler<ShapeFunction, IntegrationMethod,
DisplacementDim>(
- e, local_matrix_size, integration_order, process_data)
+ e, local_matrix_size, is_axially_symmetric, integration_order,
+ process_data)
{
}
};
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index 4a7f13ab5a838e67814e00d110a878b797f9ac6b..27bf988a1eb6e9241b63cefefaf969b03aac5edc 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -62,7 +62,8 @@ private:
ProcessLib::SmallDeformation::createLocalAssemblers<DisplacementDim,
LocalAssemblerData>(
mesh.getDimension(), mesh.getElements(), dof_table,
- integration_order, _local_assemblers, _process_data);
+ _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+ _process_data);
// TODO move the two data members somewhere else.
// for extrapolation of secondary variables
@@ -88,11 +89,31 @@ private:
getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYY));
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_zz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtSigmaZZ));
+
Base::_secondary_variables.addSecondaryVariable(
"sigma_xy", 1,
makeExtrapolator(
getExtrapolator(), _local_assemblers,
- &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY));
+ &SmallDeformationLocalAssemblerInterface::getIntPtSigmaZZ));
+
+ if (DisplacementDim == 3) {
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtSigmaYZ));
+ }
}
void assembleConcreteProcess(const double t, GlobalVector const& x,
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 86bc01e390b78a29a14b073e0afdebe749d0ce6d..43597ba518f68c323c1cb6e4f120d71619199235 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -110,12 +110,13 @@ template <typename ShapeFunction_, typename IntegrationMethod_,
TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
TESLocalAssembler(MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
unsigned const integration_order,
AssemblyParams const& asm_params)
: _integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod_, GlobalDim>(
- e, _integration_method)),
+ e, is_axially_symmetric, _integration_method)),
_d(asm_params,
// TODO narrowing conversion
static_cast<const unsigned>(
@@ -154,8 +155,8 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::assemble(
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const weight = wp.getWeight();
- _d.assembleIntegrationPoint(ip, local_x, sm.N, sm.dNdx, sm.J, sm.detJ,
- weight, local_M, local_K, local_b);
+ _d.assembleIntegrationPoint(ip, local_x, sm, weight, local_M, local_K,
+ local_b);
}
if (_d.getAssemblyParameters().output_element_matrices)
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index 86adf8f649e0607fff1f71011761070c4236ab24..5b9b4cd232d7230dff36f46680f4a8c824437ab5 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -66,6 +66,7 @@ public:
TESLocalAssembler(MeshLib::Element const& e,
std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
unsigned const integration_order,
AssemblyParams const& asm_params);
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl.h b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
index 26d529132dd90974be9167879366a43072799e92..0e5c5705de04ca11446c3dd923e6511d73d487e4 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner-impl.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
@@ -191,10 +191,7 @@ template <typename Traits>
void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
const unsigned int_pt,
const std::vector<double>& localX,
- typename Traits::ShapeMatrices::ShapeType const& smN,
- typename Traits::ShapeMatrices::DxShapeType const& /*smDNdx*/,
- typename Traits::ShapeMatrices::JacobianType const& /*smJ*/,
- const double /*smDetJ*/)
+ typename Traits::ShapeMatrices const& sm)
{
#ifndef NDEBUG
// fill local data with garbage to aid in debugging
@@ -202,7 +199,7 @@ void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
std::numeric_limits<double>::quiet_NaN();
#endif
- NumLib::shapeFunctionInterpolate(localX, smN, _d.p, _d.T,
+ NumLib::shapeFunctionInterpolate(localX, sm.N, _d.p, _d.T,
_d.vapour_mass_fraction);
// pre-compute certain properties
@@ -222,20 +219,16 @@ template <typename Traits>
void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
unsigned integration_point,
std::vector<double> const& localX,
- const typename Traits::ShapeMatrices::ShapeType& smN,
- const typename Traits::ShapeMatrices::DxShapeType& smDNdx,
- const typename Traits::ShapeMatrices::JacobianType& smJ,
- const double smDetJ,
+ typename Traits::ShapeMatrices const& sm,
const double weight,
Eigen::Map<typename Traits::LocalMatrix>& local_M,
Eigen::Map<typename Traits::LocalMatrix>& local_K,
Eigen::Map<typename Traits::LocalVector>& local_b)
{
- preEachAssembleIntegrationPoint(integration_point, localX, smN, smDNdx, smJ,
- smDetJ);
+ preEachAssembleIntegrationPoint(integration_point, localX, sm);
- auto const N = smDNdx.cols(); // number of integration points
- auto const D = smDNdx.rows(); // global dimension: 1, 2 or 3
+ auto const N = sm.dNdx.cols(); // number of integration points
+ auto const D = sm.dNdx.rows(); // global dimension: 1, 2 or 3
assert(N * NODAL_DOF == local_M.cols());
@@ -246,12 +239,12 @@ void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
auto const contentCoeffMat = getContentCoeffMatrix(integration_point);
// calculate velocity
- assert((unsigned)smDNdx.rows() == _d.velocity.size() &&
- (unsigned)smDNdx.cols() == _d.velocity[0].size());
+ assert((unsigned)sm.dNdx.rows() == _d.velocity.size() &&
+ (unsigned)sm.dNdx.cols() == _d.velocity[0].size());
auto const velocity =
(Traits::blockDimDim(laplaceCoeffMat, 0, 0, D, D) *
- (smDNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
+ (sm.dNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
N) // grad_p
/
-_d.rho_GR))
@@ -263,26 +256,27 @@ void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
_d.velocity[d][integration_point] = velocity[d];
}
- auto const detJ_w_NT = (smDetJ * weight * smN.transpose()).eval();
- auto const detJ_w_NT_N = (detJ_w_NT * smN).eval();
- assert(detJ_w_NT_N.rows() == N && detJ_w_NT_N.cols() == N);
+ auto const detJ_w_im_NT =
+ (sm.detJ * weight * sm.integralMeasure * sm.N.transpose()).eval();
+ auto const detJ_w_im_NT_N = (detJ_w_im_NT * sm.N).eval();
+ assert(detJ_w_im_NT_N.rows() == N && detJ_w_im_NT_N.cols() == N);
- auto const detJ_w_NT_vT_dNdx =
- (detJ_w_NT * velocity.transpose() * smDNdx).eval();
- assert(detJ_w_NT_vT_dNdx.rows() == N && detJ_w_NT_vT_dNdx.cols() == N);
+ auto const detJ_w_im_NT_vT_dNdx =
+ (detJ_w_im_NT * velocity.transpose() * sm.dNdx).eval();
+ assert(detJ_w_im_NT_vT_dNdx.rows() == N && detJ_w_im_NT_vT_dNdx.cols() == N);
for (unsigned r = 0; r < NODAL_DOF; ++r)
{
for (unsigned c = 0; c < NODAL_DOF; ++c)
{
Traits::blockShpShp(local_K, N * r, N * c, N, N).noalias() +=
- smDetJ * weight * smDNdx.transpose() *
+ sm.detJ * weight * sm.integralMeasure * sm.dNdx.transpose() *
Traits::blockDimDim(laplaceCoeffMat, D * r, D * c, D, D) *
- smDNdx // end Laplacian part
- + detJ_w_NT_N * contentCoeffMat(r, c) +
- detJ_w_NT_vT_dNdx * advCoeffMat(r, c);
+ sm.dNdx // end Laplacian part
+ + detJ_w_im_NT_N * contentCoeffMat(r, c) +
+ detJ_w_im_NT_vT_dNdx * advCoeffMat(r, c);
Traits::blockShpShp(local_M, N * r, N * c, N, N).noalias() +=
- detJ_w_NT_N * massCoeffMat(r, c);
+ detJ_w_im_NT_N * massCoeffMat(r, c);
}
}
@@ -291,7 +285,8 @@ void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
for (unsigned r = 0; r < NODAL_DOF; ++r)
{
Traits::blockShp(local_b, N * r, N).noalias() +=
- rhsCoeffVector(r) * smN.transpose() * smDetJ * weight;
+ rhsCoeffVector(r) * sm.N.transpose() * sm.detJ * weight *
+ sm.integralMeasure;
}
}
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.h b/ProcessLib/TES/TESLocalAssemblerInner.h
index d1bffda68f15f7ef0f525e8027281b58628a2e4d..cdb6f9ec942c2d36dce11c39c37704a986d09d25 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner.h
@@ -34,10 +34,7 @@ public:
void assembleIntegrationPoint(
unsigned integration_point,
std::vector<double> const& localX,
- typename Traits::ShapeMatrices::ShapeType const& smN,
- typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
- typename Traits::ShapeMatrices::JacobianType const& smJ,
- const double smDetJ,
+ typename Traits::ShapeMatrices const& sm,
const double weight,
Eigen::Map<typename Traits::LocalMatrix>& local_M,
Eigen::Map<typename Traits::LocalMatrix>& local_K,
@@ -64,10 +61,7 @@ private:
void preEachAssembleIntegrationPoint(
const unsigned int_pt,
std::vector<double> const& localX,
- typename Traits::ShapeMatrices::ShapeType const& smN,
- typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
- typename Traits::ShapeMatrices::JacobianType const& smJ,
- const double smDetJ);
+ typename Traits::ShapeMatrices const& sm);
void initReaction(const unsigned int_pt);
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index b8532baeb5049ff2243809077597c25753e7e749..ed754c9485b3686f126aba286001c968338dd23a 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -140,8 +140,8 @@ void TESProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh, unsigned const integration_order)
{
ProcessLib::createLocalAssemblers<TESLocalAssembler>(
- mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
- _local_assemblers, _assembly_params);
+ mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _assembly_params);
initializeSecondaryVariables();
}
diff --git a/ProcessLib/Utils/CreateLocalAssemblers.h b/ProcessLib/Utils/CreateLocalAssemblers.h
index e14aae36025f633d0d71cc71af825c0a23d45408..5e96d2f32b2260c21d2b264ac097b00b56dfcdbb 100644
--- a/ProcessLib/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/Utils/CreateLocalAssemblers.h
@@ -31,7 +31,6 @@ template<unsigned GlobalDim,
void createLocalAssemblers(
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<MeshLib::Element*> const& mesh_elements,
- unsigned const integration_order,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args
)
@@ -51,11 +50,8 @@ void createLocalAssemblers(
DBUG("Calling local assembler builder for all mesh elements.");
GlobalExecutor::transformDereferenced(
- initializer,
- mesh_elements,
- local_assemblers,
- integration_order,
- std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ initializer, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
}
} // namespace detail
@@ -80,7 +76,6 @@ void createLocalAssemblers(
const unsigned dimension,
std::vector<MeshLib::Element*> const& mesh_elements,
NumLib::LocalToGlobalIndexMap const& dof_table,
- unsigned const integration_order,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
ExtraCtorArgs&&... extra_ctor_args
)
@@ -90,25 +85,19 @@ void createLocalAssemblers(
switch (dimension)
{
case 1:
- detail::createLocalAssemblers<
- 1, LocalAssemblerImplementation>(
- dof_table, mesh_elements, integration_order,
- local_assemblers,
- std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ detail::createLocalAssemblers<1, LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
case 2:
- detail::createLocalAssemblers<
- 2, LocalAssemblerImplementation>(
- dof_table, mesh_elements, integration_order,
- local_assemblers,
- std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ detail::createLocalAssemblers<2, LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
case 3:
- detail::createLocalAssemblers<
- 3, LocalAssemblerImplementation>(
- dof_table, mesh_elements, integration_order,
- local_assemblers,
- std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ detail::createLocalAssemblers<3, LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
break;
default:
OGS_FATAL("Meshes with dimension greater than three are not supported.");
diff --git a/ProcessLib/Utils/InitShapeMatrices.h b/ProcessLib/Utils/InitShapeMatrices.h
index 6fbf385775f96d6e54298ff401d3a926b6ec62c7..78e3b0459dad0633f51f9f1b8c52719fe2946f73 100644
--- a/ProcessLib/Utils/InitShapeMatrices.h
+++ b/ProcessLib/Utils/InitShapeMatrices.h
@@ -20,7 +20,8 @@ namespace ProcessLib
template <typename ShapeFunction, typename ShapeMatricesType,
typename IntegrationMethod, unsigned GlobalDim>
std::vector<typename ShapeMatricesType::ShapeMatrices> initShapeMatrices(
- MeshLib::Element const& e, IntegrationMethod const& integration_method)
+ MeshLib::Element const& e, bool is_axially_symmetric,
+ IntegrationMethod const& integration_method)
{
std::vector<typename ShapeMatricesType::ShapeMatrices> shape_matrices;
@@ -37,12 +38,25 @@ std::vector<typename ShapeMatricesType::ShapeMatrices> initShapeMatrices(
ShapeFunction::NPOINTS);
fe.computeShapeFunctions(
integration_method.getWeightedPoint(ip).getCoords(),
- shape_matrices[ip], GlobalDim);
+ shape_matrices[ip], GlobalDim, is_axially_symmetric);
}
return shape_matrices;
}
+template <typename ShapeFunction, typename ShapeMatricesType>
+double interpolateXCoordinate(
+ MeshLib::Element const& e,
+ typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
+{
+ using FemType = NumLib::TemplateIsoparametric<
+ ShapeFunction, ShapeMatricesType>;
+
+ FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
+
+ return fe.interpolateZerothCoordinate(N);
+}
+
} // ProcessLib
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
index d1818cbcf0d52075b04e5e6c1d6fb43989259e16..fe7bdfa2bbdf51acb1e8bc9176a1f0b0a29ba130 100644
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ b/ProcessLib/Utils/LocalDataInitializer.h
@@ -249,7 +249,6 @@ public:
std::size_t const id,
MeshLib::Element const& mesh_item,
LADataIntfPtr& data_ptr,
- unsigned const integration_order,
ConstructorArgs&&... args) const
{
auto const type_idx = std::type_index(typeid(mesh_item));
@@ -257,9 +256,8 @@ public:
if (it != _builder.end()) {
auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
- data_ptr = it->second(
- mesh_item, num_local_dof, integration_order,
- std::forward<ConstructorArgs>(args)...);
+ data_ptr = it->second(mesh_item, num_local_dof,
+ std::forward<ConstructorArgs>(args)...);
} else {
OGS_FATAL("You are trying to build a local assembler for an unknown mesh element type (%s)."
" Maybe you have disabled this mesh element type in your build configuration.",
@@ -271,7 +269,6 @@ private:
using LADataBuilder = std::function<LADataIntfPtr(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
- unsigned const integration_order,
ConstructorArgs&&...
)>;
@@ -291,12 +288,11 @@ private:
{
return [](MeshLib::Element const& e,
std::size_t const local_matrix_size,
- unsigned const integration_order,
ConstructorArgs&&... args)
{
return LADataIntfPtr{
new LAData<ShapeFct>{
- e, local_matrix_size, integration_order,
+ e, local_matrix_size,
std::forward<ConstructorArgs>(args)...
}};
};
diff --git a/Tests/Data b/Tests/Data
index 510f6ba3f243b48ef4e170c86b69dc235dae4bff..6f8577bb599780d0e4fc99580da282e928173dc6 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 510f6ba3f243b48ef4e170c86b69dc235dae4bff
+Subproject commit 6f8577bb599780d0e4fc99580da282e928173dc6
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index aee0e030dd1ce15a3efc358366950b17189babb7..20f9c9a70a8dc04e81129ec357d4d4d2ebd0e24c 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -73,11 +73,13 @@ class LocalAssemblerData : public LocalAssemblerDataInterface
public:
LocalAssemblerData(MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
unsigned const integration_order)
: _shape_matrices(
ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
- e, IntegrationMethod{integration_order})),
+ e, is_axially_symmetric,
+ IntegrationMethod{integration_order})),
_int_pt_values(_shape_matrices.size())
{
}
@@ -142,7 +144,10 @@ public:
// and the variable has exactly one component.
_extrapolator.reset(new ExtrapolatorImplementation(*_dof_table));
- createAssemblers(mesh);
+ // createAssemblers(mesh);
+ ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+ mesh.getDimension(), mesh.getElements(), *_dof_table,
+ _local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
}
void interpolateNodalValuesToIntegrationPoints(
@@ -174,36 +179,6 @@ public:
}
private:
- void createAssemblers(MeshLib::Mesh const& mesh)
- {
- switch (mesh.getDimension())
- {
- case 1: createLocalAssemblers<1>(mesh); break;
- case 2: createLocalAssemblers<2>(mesh); break;
- case 3: createLocalAssemblers<3>(mesh); break;
- default: assert(false);
- }
- }
-
- template<unsigned GlobalDim>
- void createLocalAssemblers(MeshLib::Mesh const& mesh)
- {
- using LocalDataInitializer = ProcessLib::LocalDataInitializer<
- LocalAssembler, LocalAssemblerData,
- GlobalDim>;
-
- _local_assemblers.resize(mesh.getNumberOfElements());
-
- LocalDataInitializer initializer(*_dof_table);
-
- DBUG("Calling local assembler builder for all mesh elements.");
- GlobalExecutor::transformDereferenced(
- initializer,
- mesh.getElements(),
- _local_assemblers,
- _integration_order);
- }
-
unsigned const _integration_order;
MeshLib::MeshSubset _mesh_subset_all_nodes;
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 6b031f20d811ea2d7af66e85c115c87e09a97795..b03cecb8107aa4f38b138fe71d47b5a1c181a879 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -195,7 +195,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
fe.template computeShapeFunctions<ShapeMatrixType::N_J>(
- wp.getCoords(), shape, this->global_dim);
+ wp.getCoords(), shape, this->global_dim, false);
M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M;
@@ -220,7 +220,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(
- wp.getCoords(), shape, this->global_dim);
+ wp.getCoords(), shape, this->global_dim, false);
K.noalias() += shape.dNdx.transpose() * this->globalD * shape.dNdx * shape.detJ * wp.getWeight();
}
//std::cout << "K=\n" << K << std::endl;
@@ -246,7 +246,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape, this->global_dim);
+ fe.computeShapeFunctions(wp.getCoords(), shape, this->global_dim, false);
M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
K.noalias() += shape.dNdx.transpose() * (this->globalD * shape.dNdx) * shape.detJ * wp.getWeight();
}
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index 7e507407f80d9ad8dd79ca19004d47acf53ef815..ff1a1b32b1bd86838a3585972b671ae88aad3bed 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -71,7 +71,7 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
finite_element.computeShapeFunctions(
integration_method.getWeightedPoint(0).getCoords(),
- shape_matrix, ShapeFunction::DIM);
+ shape_matrix, ShapeFunction::DIM, false);
ASSERT_EQ(2, shape_matrix.N.size());
// actual test