OGS
SurfaceFluxLocalAssembler.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Geometry>
14
22
23namespace ProcessLib
24{
26{
27public:
29
30 virtual void integrate(
31 std::size_t const element_id,
32 std::vector<GlobalVector*> const& x,
34 double const t,
35 MeshLib::Mesh const& bulk_mesh,
36 std::function<Eigen::Vector3d(std::size_t const,
37 MathLib::Point3d const&, double const,
38 std::vector<GlobalVector*> const&)>) = 0;
39};
40
41template <typename ShapeFunction, int GlobalDim>
44{
45protected:
49
50public:
62 MeshLib::Element const& surface_element,
63 std::size_t /*const local_matrix_size*/,
64 NumLib::GenericIntegrationMethod const& integration_method,
65 bool const is_axially_symmetric,
66 MeshLib::PropertyVector<std::size_t> const& bulk_element_ids,
67 MeshLib::PropertyVector<std::size_t> const& bulk_face_ids)
68 : _surface_element(surface_element),
69 _integration_method(integration_method),
70 _bulk_element_id(bulk_element_ids[surface_element.getID()]),
71 _bulk_face_id(bulk_face_ids[surface_element.getID()])
72 {
73 auto const n_integration_points =
75
76 auto const shape_matrices =
79 _surface_element, is_axially_symmetric, _integration_method);
80 for (std::size_t ip = 0; ip < n_integration_points; ++ip)
81 {
83 shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure);
84 }
85 }
86
100 void integrate(std::size_t const element_id,
101 std::vector<GlobalVector*> const& x,
102 MeshLib::PropertyVector<double>& specific_flux,
103 double const t,
104 MeshLib::Mesh const& bulk_mesh,
105 std::function<Eigen::Vector3d(
106 std::size_t const, MathLib::Point3d const&, double const,
107 std::vector<GlobalVector*> const&)>
108 getFlux) override
109 {
110 auto const& bulk_element = *bulk_mesh.getElement(_bulk_element_id);
111 auto const surface_element_normal =
112 getSurfaceNormal(_surface_element, bulk_element);
113
114 double element_area = 0.0;
115 std::size_t const n_integration_points =
117 // specific_flux[id_of_element] +=
118 // int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma /
119 // element_area
120 for (std::size_t ip(0); ip < n_integration_points; ip++)
121 {
122 auto const& wp = _integration_method.getWeightedPoint(ip);
123
124 auto const bulk_element_point =
126 auto const bulk_flux =
127 getFlux(_bulk_element_id, bulk_element_point, t, x);
128 for (int component_id(0);
129 component_id < specific_flux.getNumberOfGlobalComponents();
130 ++component_id)
131 {
132 // TODO find solution for 2d case
133 double const bulk_grad_times_normal(
134 Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
135 bulk_flux.size())
136 .dot(surface_element_normal));
137
138 specific_flux.getComponent(element_id, component_id) +=
139 bulk_grad_times_normal * _detJ_times_integralMeasure[ip] *
140 wp.getWeight();
141 }
142 element_area += _detJ_times_integralMeasure[ip] * wp.getWeight();
143 }
144 for (int component_id(0);
145 component_id < specific_flux.getNumberOfGlobalComponents();
146 ++component_id)
147 {
148 specific_flux.getComponent(element_id, component_id) /=
149 element_area;
150 }
151 }
152
153private:
154 static Eigen::Vector3d getSurfaceNormal(
155 MeshLib::Element const& surface_element,
156 MeshLib::Element const& bulk_element)
157 {
158 Eigen::Vector3d surface_element_normal;
159
160 if (surface_element.getGeomType() == MeshLib::MeshElemType::LINE)
161 {
162 auto const bulk_normal =
164 auto const& l0 = surface_element.getNode(0)->asEigenVector3d();
165 auto const& l1 = surface_element.getNode(1)->asEigenVector3d();
166 Eigen::Vector3d const line = l1 - l0;
167 surface_element_normal = line.cross(bulk_normal);
168 }
169 else
170 {
171 surface_element_normal =
173 }
174 surface_element_normal.normalize();
175 // At the moment (2016-09-28) the surface normal is not oriented
176 // according to the right hand rule. Thus for an intuitive flux
177 // output, i.e., inflow has positive sign, outflow has negative
178 // sign, the normal must not be multiplied by -1.
179 return surface_element_normal;
180 };
181
183
184 std::vector<double> _detJ_times_integralMeasure;
185
187 std::size_t const _bulk_element_id;
188 std::size_t const _bulk_face_id;
189};
190
191} // namespace ProcessLib
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
virtual MeshElemType getGeomType() const =0
virtual const Node * getNode(unsigned idx) const =0
static Eigen::Vector3d getSurfaceNormal(Element const &e)
Returns the surface normal of a 2D element.
Definition FaceRule.cpp:40
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:94
int getNumberOfGlobalComponents() const
PROP_VAL_TYPE & getComponent(std::size_t tuple_index, int component)
Returns the value for the given component stored in the given tuple.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
virtual void integrate(std::size_t const element_id, std::vector< GlobalVector * > const &x, MeshLib::PropertyVector< double > &specific_flux, double const t, MeshLib::Mesh const &bulk_mesh, std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)>)=0
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
static Eigen::Vector3d getSurfaceNormal(MeshLib::Element const &surface_element, MeshLib::Element const &bulk_element)
void integrate(std::size_t const element_id, std::vector< GlobalVector * > const &x, MeshLib::PropertyVector< double > &specific_flux, double const t, MeshLib::Mesh const &bulk_mesh, std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)> getFlux) override
NumLib::GenericIntegrationMethod const & _integration_method
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
SurfaceFluxLocalAssembler(MeshLib::Element const &surface_element, std::size_t, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, MeshLib::PropertyVector< std::size_t > const &bulk_element_ids, MeshLib::PropertyVector< std::size_t > const &bulk_face_ids)
typename ShapeMatricesType::NodalVectorType NodalVectorType
MathLib::Point3d getBulkElementPoint(MeshLib::CellType const bulk_element_cell_type, std::size_t const bulk_face_id, MathLib::WeightedPoint const &point_on_face)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
@ N_J
calculates N, dNdr, J, and detJ
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< ShapeFunction::NPOINTS > NodalVectorType