https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: f64b0e204baab829f407ca0451b99e23749b8cc2 authored by Lars Bilke on 15 March 2021, 09:36:53 UTC
Bump Boost to 1.69.0 for Clang compatibility.
Bump Boost to 1.69.0 for Clang compatibility.
Tip revision: f64b0e2
SurfaceFluxLocalAssembler.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include "MeshLib/Elements/Elements.h"
#include "MeshLib/Elements/FaceRule.h"
#include "MeshLib/Elements/MapBulkElementPoint.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Process.h"
namespace ProcessLib
{
class SurfaceFluxLocalAssemblerInterface
{
public:
virtual ~SurfaceFluxLocalAssemblerInterface() = default;
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;
};
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
class SurfaceFluxLocalAssembler final
: public SurfaceFluxLocalAssemblerInterface
{
protected:
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
public:
/// Precomputes the shape matrices for a given surface element.
/// @param surface_element The surface element used for precomputing the
/// shape matrices used later on for the integration.
/// @param is_axially_symmetric corrects integration measure for cylinder
/// coordinates.
/// @param bulk_element_ids The id of the corresponding element in the bulk
/// mesh.
/// @param bulk_face_ids The id of the corresponding face in the bulk
/// element.
/// @param integration_order the order of the integration
SurfaceFluxLocalAssembler(
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)
: _surface_element(surface_element),
_integration_method(integration_order),
_bulk_element_id(bulk_element_ids[surface_element.getID()]),
_bulk_face_id(bulk_face_ids[surface_element.getID()])
{
auto const n_integration_points =
_integration_method.getNumberOfPoints();
auto const shape_matrices =
NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
GlobalDim, NumLib::ShapeMatrixType::N_J>(
_surface_element, is_axially_symmetric, _integration_method);
for (std::size_t ip = 0; ip < n_integration_points; ++ip)
{
_detJ_times_integralMeasure.push_back(
shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure);
}
}
/// Integration for the element with the id \c element_id.
/// @param element_id id of the element
/// @param x The global vector containing the values for numerical
/// integration.
/// @param specific_flux PropertyVector the integration result will be
/// stored into, where specific_flux has the same number of entries like
/// mesh elements exists.
/// @param t The integration is performed at the time t.
/// @param bulk_mesh Stores a reference to the bulk mesh that is needed to
/// fetch the information for the integration over the surface mesh.
/// @param getFlux function that calculates the flux in the integration
/// points of the face elements of the bulk element that belongs to the
/// surface.
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
{
auto get_surface_normal =
[this, &bulk_mesh](
MeshLib::Element const& surface_element) -> Eigen::Vector3d {
Eigen::Vector3d surface_element_normal;
if (surface_element.getGeomType() == MeshLib::MeshElemType::LINE)
{
auto const bulk_normal = MeshLib::FaceRule::getSurfaceNormal(
bulk_mesh.getElements()[_bulk_element_id]);
auto const l0 = Eigen::Map<Eigen::Vector3d const>(
_surface_element.getNode(0)->getCoords());
auto const l1 = Eigen::Map<Eigen::Vector3d const>(
_surface_element.getNode(1)->getCoords());
Eigen::Vector3d const line = l1 - l0;
surface_element_normal = line.cross(bulk_normal);
}
else
{
surface_element_normal =
MeshLib::FaceRule::getSurfaceNormal(&surface_element);
}
surface_element_normal.normalize();
// At the moment (2016-09-28) the surface normal is not oriented
// according to the right hand rule. Thus for an intuitive flux
// output, i.e., inflow has positive sign, outflow has negative
// sign, the normal must not be multiplied by -1.
return surface_element_normal;
};
auto const surface_element_normal =
get_surface_normal(_surface_element);
double element_area = 0.0;
std::size_t const n_integration_points =
_integration_method.getNumberOfPoints();
// specific_flux[id_of_element] +=
// int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma /
// element_area
for (std::size_t ip(0); ip < n_integration_points; ip++)
{
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const bulk_element_point = MeshLib::getBulkElementPoint(
bulk_mesh, _bulk_element_id, _bulk_face_id, wp);
auto const bulk_flux =
getFlux(_bulk_element_id, bulk_element_point, t, x);
for (int component_id(0);
component_id < specific_flux.getNumberOfGlobalComponents();
++component_id)
{
// TODO find solution for 2d case
double const bulk_grad_times_normal(
Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
bulk_flux.size())
.dot(surface_element_normal));
specific_flux.getComponent(element_id, component_id) +=
bulk_grad_times_normal * _detJ_times_integralMeasure[ip] *
wp.getWeight();
}
element_area += _detJ_times_integralMeasure[ip] * wp.getWeight();
}
for (int component_id(0);
component_id < specific_flux.getNumberOfGlobalComponents();
++component_id)
{
specific_flux.getComponent(element_id, component_id) /=
element_area;
}
}
private:
MeshLib::Element const& _surface_element;
std::vector<double> _detJ_times_integralMeasure;
IntegrationMethod const _integration_method;
std::size_t _bulk_element_id;
std::size_t _bulk_face_id;
};
} // namespace ProcessLib