From 9d8a9cd4457cb78601b7982a7ca429e55d2d757c Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Tue, 26 Jun 2018 09:16:39 +0200
Subject: [PATCH] [PL/ConstraintDirichletBC] Mv calc. of normal to constr.
---
.../ConstraintDirichletBoundaryCondition.cpp | 2 +-
...DirichletBoundaryConditionLocalAssembler.h | 52 +++++++++++--------
2 files changed, 31 insertions(+), 23 deletions(-)
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
index b85cd94d843..7f74f4a7662 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -116,7 +116,7 @@ ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
ConstraintDirichletBoundaryConditionLocalAssembler>(
_bulk_mesh.getDimension(), _bc_mesh.getElements(), *_dof_table_boundary,
shape_function_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
- _integration_order, _bulk_ids);
+ _integration_order, bulk_mesh, _bulk_ids);
}
void ConstraintDirichletBoundaryCondition::preTimestep(double t,
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
index d90e28a9285..dee62b81603 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -23,6 +23,7 @@
namespace ProcessLib
{
+
struct IntegrationPointData final
{
IntegrationPointData(double const& detJ,
@@ -70,11 +71,13 @@ public:
/// @param is_axially_symmetric Corrects integration measure for cylinder
/// coordinates.
/// @param integration_order The order of the integration.
+ /// @param bulk_mesh The bulk mesh the process is defined on.
/// @param bulk_ids Pairs of bulk element ids and bulk element face ids.
ConstraintDirichletBoundaryConditionLocalAssembler(
MeshLib::Element const& surface_element,
std::size_t const local_matrix_size,
bool const is_axially_symmetric, unsigned const integration_order,
+ MeshLib::Mesh const& bulk_mesh,
std::vector<std::pair<std::size_t, unsigned>> bulk_ids)
: _surface_element(surface_element),
_integration_method(integration_order),
@@ -82,6 +85,30 @@ public:
_bulk_face_id(bulk_ids[_surface_element.getID()].second)
{
(void)local_matrix_size; // unused, but needed for the interface
+
+ if (_surface_element.getDimension() < 2)
+ {
+ auto const bulk_element_normal =
+ MeshLib::FaceRule::getSurfaceNormal(
+ bulk_mesh.getElement(_bulk_element_id));
+ MathLib::Vector3 const edge_vector(*_surface_element.getNode(0),
+ *_surface_element.getNode(1));
+ _surface_element_normal =
+ MathLib::crossProduct(bulk_element_normal, edge_vector);
+ }
+ else
+ {
+ _surface_element_normal =
+ MeshLib::FaceRule::getSurfaceNormal(&_surface_element);
+ }
+
+ _surface_element_normal.normalize();
+ // At the moment (2018-04-26) the surface normal is not oriented
+ // according to the right hand rule
+ // for correct results it is necessary to multiply the normal with
+ // -1
+ _surface_element_normal *= -1;
+
using FemType =
NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
@@ -104,6 +131,7 @@ public:
fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N_J>(
_integration_method.getWeightedPoint(ip).getCoords(),
shape_matrices[ip], GlobalDim, is_axially_symmetric);
+
auto const& wp = _integration_method.getWeightedPoint(ip);
_ip_data.emplace_back(shape_matrices[ip].detJ,
shape_matrices[ip].integralMeasure,
@@ -124,27 +152,6 @@ public:
std::size_t const, MathLib::Point3d const&, double const,
GlobalVector const&)> const& getFlux) override
{
- MathLib::Vector3 surface_element_normal;
- if (_surface_element.getDimension() < 2)
- {
- auto const bulk_element_normal =
- MeshLib::FaceRule::getSurfaceNormal(
- bulk_mesh.getElement(_bulk_element_id));
- MathLib::Vector3 const edge_vector(*_surface_element.getNode(0),
- *_surface_element.getNode(1));
- surface_element_normal =
- MathLib::crossProduct(bulk_element_normal, edge_vector);
- } else {
- surface_element_normal =
- MeshLib::FaceRule::getSurfaceNormal(&_surface_element);
- }
-
- surface_element_normal.normalize();
- // At the moment (2018-04-26) the surface normal is not oriented
- // according to the right hand rule
- // for correct results it is necessary to multiply the normal with -1
- surface_element_normal *= -1;
-
auto const n_integration_points =
_integration_method.getNumberOfPoints();
// integrated_value +=
@@ -163,7 +170,7 @@ public:
Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
bulk_flux.size())
.dot(Eigen::Map<Eigen::RowVectorXd const>(
- surface_element_normal.getCoords(), 3)));
+ _surface_element_normal.getCoords(), 3)));
integrated_value +=
bulk_grad_times_normal *
@@ -180,6 +187,7 @@ private:
IntegrationMethod const _integration_method;
std::size_t const _bulk_element_id;
unsigned const _bulk_face_id;
+ MathLib::Vector3 _surface_element_normal;
};
} // ProcessLib
--
GitLab