diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 263e962cd0d990c864ecd87a0e13aa7e9af886ff..0f09c01d8a8d7d5be1a4ee471ef20da288defdfd 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -54,28 +54,14 @@ public:
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();
- SpatialPosition pos;
- pos.setElementID(id);
-
- using FemType =
- NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
- FemType fe(
- static_cast<const typename ShapeFunction::MeshElement&>(Base::_element));
-
- // Get element nodes for the interpolation from nodes to
- // integration point.
+ // Get element nodes for the interpolation from nodes to integration
+ // point.
NodalVectorType parameter_node_values =
_neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t);
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- pos.setIntegrationPoint(ip);
auto const& ip_data = Base::_ns_and_weights[ip];
- auto const& wp = Base::_integration_method.getWeightedPoint(ip);
-
- MathLib::TemplatePoint<double, 3> coords_ip(
- fe.interpolateCoordinates(ip_data.N));
- pos.setCoordinates(coords_ip);
_local_rhs.noalias() += ip_data.N *
parameter_node_values.dot(ip_data.N) *
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
index 7f1315dced38b34461fd55ae631de1b4d8e2687d..368a45ef2e88a1f4ac8025411eb383f054109a61 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
@@ -21,21 +21,19 @@ namespace ProcessLib
{
namespace NormalTractionBoundaryCondition
{
-template <typename ShapeMatricesTypeDisplacement, int GlobalDim, int NPoints>
+template <typename ShapeMatricesType>
struct IntegrationPointData final
{
IntegrationPointData(
- typename ShapeMatricesTypeDisplacement::template VectorType<
- NPoints * GlobalDim> const& Nu_times_n_,
+ typename ShapeMatricesType::ShapeMatrices::ShapeType const N_,
+ typename ShapeMatricesType::GlobalDimVectorType const n_,
double const integration_weight_)
- : Nu_times_n(Nu_times_n_), integration_weight(integration_weight_)
+ : N(N_), n(n_), integration_weight(integration_weight_)
{
}
- // Shape matrix (for displacement) times element's normal vector.
- typename ShapeMatricesTypeDisplacement::template VectorType<
- NPoints * GlobalDim> const Nu_times_n;
-
+ typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
+ typename ShapeMatricesType::GlobalDimVectorType const n;
double const integration_weight;
};
@@ -56,11 +54,10 @@ class NormalTractionBoundaryConditionLocalAssembler final
: public NormalTractionBoundaryConditionLocalAssemblerInterface
{
public:
- using ShapeMatricesTypeDisplacement =
+ using ShapeMatricesType =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
- using GlobalDimVectorType =
- typename ShapeMatrixPolicyType<ShapeFunctionDisplacement,
- GlobalDim>::GlobalDimVectorType;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
NormalTractionBoundaryConditionLocalAssembler(
MeshLib::Element const& e,
@@ -68,7 +65,9 @@ public:
bool const is_axially_symmetric,
unsigned const integration_order,
Parameter<double> const& pressure)
- : _integration_method(integration_order), _pressure(pressure)
+ : _integration_method(integration_order),
+ _pressure(pressure),
+ _element(e)
{
_local_rhs.setZero(local_matrix_size);
@@ -78,10 +77,9 @@ public:
_ip_data.reserve(n_integration_points);
auto const shape_matrices_u =
- initShapeMatrices<ShapeFunctionDisplacement,
- ShapeMatricesTypeDisplacement, IntegrationMethod,
- GlobalDim>(e, is_axially_symmetric,
- _integration_method);
+ initShapeMatrices<ShapeFunctionDisplacement, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ e, is_axially_symmetric, _integration_method);
GlobalDimVectorType element_normal(GlobalDim);
@@ -104,24 +102,13 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- typename ShapeMatricesTypeDisplacement::template MatrixType<
- GlobalDim, displacement_size>
- N_u = ShapeMatricesTypeDisplacement::template MatrixType<
- GlobalDim, displacement_size>::Zero(GlobalDim,
- displacement_size);
- for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
- {
- N_u.template block<1, displacement_size / GlobalDim>(
- i, i * displacement_size / GlobalDim)
- .noalias() = shape_matrices_u[ip].N;
- }
double const integration_weight =
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices_u[ip].integralMeasure *
shape_matrices_u[ip].detJ;
- _ip_data.emplace_back(N_u.transpose() * element_normal,
+ _ip_data.emplace_back(shape_matrices_u[ip].N, element_normal,
integration_weight);
}
}
@@ -137,18 +124,28 @@ public:
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
- SpatialPosition pos;
- pos.setElementID(id);
+ NodalVectorType pressure =
+ _pressure.getNodalValuesOnElement(_element, t);
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- pos.setIntegrationPoint(ip);
-
auto const& w = _ip_data[ip].integration_weight;
- auto const& Nu_times_n = _ip_data[ip].Nu_times_n;
+ auto const& N = _ip_data[ip].N;
+ auto const& n = _ip_data[ip].n;
- _local_rhs.noalias() -=
- Nu_times_n.transpose() * _pressure(t, pos)[0] * w;
+ typename ShapeMatricesType::template MatrixType<GlobalDim,
+ displacement_size>
+ N_u = ShapeMatricesType::template MatrixType<
+ GlobalDim, displacement_size>::Zero(GlobalDim,
+ displacement_size);
+ for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
+ {
+ N_u.template block<1, displacement_size / GlobalDim>(
+ i, i * displacement_size / GlobalDim)
+ .noalias() = N;
+ }
+
+ _local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w;
}
auto const indices = NumLib::getIndices(id, dof_table_boundary);
@@ -161,17 +158,16 @@ private:
static const int displacement_size =
ShapeFunctionDisplacement::NPOINTS * GlobalDim;
- std::vector<IntegrationPointData<ShapeMatricesTypeDisplacement, GlobalDim,
- ShapeFunctionDisplacement::NPOINTS>,
- Eigen::aligned_allocator<IntegrationPointData<
- ShapeMatricesTypeDisplacement, GlobalDim,
- ShapeFunctionDisplacement::NPOINTS>>>
+ std::vector<
+ IntegrationPointData<ShapeMatricesType>,
+ Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>>
_ip_data;
- typename ShapeMatricesTypeDisplacement::template VectorType<
- displacement_size>
+ typename ShapeMatricesType::template VectorType<displacement_size>
_local_rhs;
+ MeshLib::Element const& _element;
+
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
index ca87bdcd3a626274cdbf7ce592c02e6c64f30457..40437fa3317ff8863344972b7846068a8a1ccee6 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
@@ -41,9 +41,7 @@ public:
bool is_axially_symmetric,
unsigned const integration_order,
PythonBoundaryConditionData const& data)
- : Base(e, is_axially_symmetric, integration_order),
- _data(data),
- _element(e)
+ : Base(e, is_axially_symmetric, integration_order), _data(data)
{
}
@@ -57,12 +55,12 @@ public:
using FemType =
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
- &_element));
+ &this->_element));
unsigned const num_integration_points =
Base::_integration_method.getNumberOfPoints();
auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
- auto const num_nodes = _element.getNumberOfNodes();
+ auto const num_nodes = Base::_element.getNumberOfNodes();
auto const num_comp_total =
_data.dof_table_bulk.getNumberOfComponents();
@@ -85,7 +83,8 @@ public:
for (unsigned element_node_id = 0; element_node_id < num_nodes;
++element_node_id)
{
- auto const* const node = _element.getNode(element_node_id);
+ auto const* const node =
+ Base::_element.getNode(element_node_id);
auto const boundary_node_id = node->getID();
auto const bulk_node_id =
bulk_node_ids_map[boundary_node_id];
@@ -119,12 +118,12 @@ public:
for (unsigned ip = 0; ip < num_integration_points; ip++)
{
- auto const& sm = Base::_shape_matrices[ip];
- auto const coords = fe.interpolateCoordinates(sm.N);
+ auto const& N = Base::_ns_and_weights[ip].N;
+ auto const& w = Base::_ns_and_weights[ip].weight;
+ auto const coords = fe.interpolateCoordinates(N);
prim_vars =
- sm.N *
- primary_variables_mat; // Assumption: all primary variables
- // have same shape functions.
+ N * primary_variables_mat; // Assumption: all primary variables
+ // have same shape functions.
auto const flag_flux_dFlux =
_data.bc_object->getFlux(t, coords, prim_vars_data);
if (!_data.bc_object->isOverriddenNatural())
@@ -143,9 +142,7 @@ public:
auto const flux = std::get<1>(flag_flux_dFlux);
auto const& dFlux = std::get<2>(flag_flux_dFlux);
- auto const& wp = Base::_integration_method.getWeightedPoint(ip);
- auto const w = sm.detJ * wp.getWeight() * sm.integralMeasure;
- local_rhs.noalias() += sm.N * (flux * w);
+ local_rhs.noalias() += N * (flux * w);
if (static_cast<int>(dFlux.size()) != num_comp_total)
{
@@ -170,7 +167,7 @@ public:
// The assignement -= takes into account the sign convention
// of 1st-order in time ODE systems in OpenGeoSys.
local_Jac.block(top, left, width, height).noalias() -=
- sm.N.transpose() * (dFlux[comp] * w) * sm.N;
+ N.transpose() * (dFlux[comp] * w) * N;
}
}
}
@@ -193,7 +190,6 @@ public:
private:
PythonBoundaryConditionData const& _data;
- MeshLib::Element const& _element;
};
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index ad5283652d017786a601de8e7b40ad4d64ed07fa..434afb14b34402374a9f42c87d7d92cb350d0779 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -54,22 +54,22 @@ public:
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();
- SpatialPosition pos;
- pos.setElementID(id);
+ typename Base::NodalVectorType const alpha =
+ _data.alpha.getNodalValuesOnElement(Base::_element, t);
+ typename Base::NodalVectorType const u_0 =
+ _data.u_0.getNodalValuesOnElement(Base::_element, t);
for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
- pos.setIntegrationPoint(ip);
auto const& ip_data = Base::_ns_and_weights[ip];
-
- double const alpha = _data.alpha(t, pos)[0];
- double const u_0 = _data.u_0(t, pos)[0];
+ auto const& N = ip_data.N;
+ auto const& w = ip_data.weight;
// flux = alpha * ( u_0 - u )
// 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() += ip_data.N * alpha * ip_data.weight;
- _local_rhs.noalias() += ip_data.N * alpha * u_0 * ip_data.weight;
+ _local_K.diagonal().noalias() += N * alpha.dot(N) * w;
+ _local_rhs.noalias() += N * alpha.dot(N) * u_0.dot(N) * w;
}
auto const indices = NumLib::getIndices(id, dof_table_boundary);
diff --git a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
index 024b3de047d264d54b6f81c56a7f544bd4f399d1..b5e0b0112f1c4c736d1796e82f0e301e1e46cc16 100644
--- a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h
@@ -74,7 +74,8 @@ public:
NodalVectorType const coefficient_mixed_variables_node_values =
_data.coefficient_mixed_variables.getNodalValuesOnElement(
Base::_element, t);
- unsigned const n_integration_points = _local_matrix_size;
+ unsigned const n_integration_points =
+ Base::_integration_method.getNumberOfPoints();
auto const indices_current_variable =
NumLib::getIndices(mesh_item_id, dof_table_boundary);