OGS
ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim > Class Template Referencefinal

Detailed Description

template<typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >

Definition at line 28 of file PythonBoundaryConditionLocalAssembler.h.

#include <PythonBoundaryConditionLocalAssembler.h>

Inheritance diagram for ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >:
[legend]

Public Member Functions

 PythonBoundaryConditionLocalAssembler (MeshLib::Element const &e, std::size_t const, bool is_axially_symmetric, unsigned const integration_order, PythonBoundaryConditionData const &data)
 
void assemble (std::size_t const boundary_element_id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const t, std::vector< GlobalVector * > const &x, int const process_id, GlobalMatrix &, GlobalVector &b, GlobalMatrix *Jac) override
 
- Public Member Functions inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >
 GenericNaturalBoundaryConditionLocalAssembler (MeshLib::Element const &e, bool is_axially_symmetric, unsigned const integration_order)
 
- Public Member Functions inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssemblerInterface
virtual ~GenericNaturalBoundaryConditionLocalAssemblerInterface ()=default
 

Private Types

using Base = GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >
 

Private Attributes

PythonBoundaryConditionData const & _data
 

Additional Inherited Members

- Protected Types inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >
using ShapeMatricesType = ShapeMatrixPolicyType< ShapeFunction, GlobalDim >
 
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType
 
using NodalVectorType = typename ShapeMatricesType::NodalVectorType
 
- Protected Attributes inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >
IntegrationMethod const _integration_method
 
std::vector< NAndWeight, Eigen::aligned_allocator< NAndWeight > > const _ns_and_weights
 
MeshLib::Element const & _element
 

Member Typedef Documentation

◆ Base

template<typename ShapeFunction , typename IntegrationMethod , int GlobalDim>
using ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::Base = GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim>
private

Definition at line 32 of file PythonBoundaryConditionLocalAssembler.h.

Constructor & Destructor Documentation

◆ PythonBoundaryConditionLocalAssembler()

template<typename ShapeFunction , typename IntegrationMethod , int GlobalDim>
ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::PythonBoundaryConditionLocalAssembler ( MeshLib::Element const &  e,
std::size_t const  ,
bool  is_axially_symmetric,
unsigned const  integration_order,
PythonBoundaryConditionData const &  data 
)
inline

Definition at line 36 of file PythonBoundaryConditionLocalAssembler.h.

42  : Base(e, is_axially_symmetric, integration_order), _data(data)
43  {
44  }
GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim > Base

Member Function Documentation

◆ assemble()

template<typename ShapeFunction , typename IntegrationMethod , int GlobalDim>
void ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::assemble ( std::size_t const  boundary_element_id,
NumLib::LocalToGlobalIndexMap const &  dof_table_boundary,
double const  t,
std::vector< GlobalVector * > const &  x,
int const  process_id,
GlobalMatrix ,
GlobalVector b,
GlobalMatrix Jac 
)
inlineoverridevirtual

Implements ProcessLib::GenericNaturalBoundaryConditionLocalAssemblerInterface.

Definition at line 46 of file PythonBoundaryConditionLocalAssembler.h.

51  {
52  using ShapeMatricesType =
55  ShapeFunction, ShapeMatricesType>(Base::_element);
56 
57  unsigned const num_integration_points =
58  Base::_integration_method.getNumberOfPoints();
59  auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
60  auto const num_nodes = Base::_element.getNumberOfNodes();
61  auto const num_comp_total =
63 
64  auto const& bulk_node_ids_map =
66  .template getPropertyVector<std::size_t>(
67  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
68 
69  // gather primary variables
70  Eigen::MatrixXd primary_variables_mat(num_nodes, num_comp_total);
71  for (int var = 0; var < num_var; ++var)
72  {
73  auto const num_comp =
75  for (int comp = 0; comp < num_comp; ++comp)
76  {
77  auto const global_component =
78  dof_table_boundary.getGlobalComponent(var, comp);
79 
80  for (unsigned element_node_id = 0; element_node_id < num_nodes;
81  ++element_node_id)
82  {
83  auto const* const node =
84  Base::_element.getNode(element_node_id);
85  auto const boundary_node_id = node->getID();
86  auto const bulk_node_id =
87  bulk_node_ids_map[boundary_node_id];
90  bulk_node_id};
91  auto const dof_idx =
92  _data.dof_table_bulk.getGlobalIndex(loc, var, comp);
93  if (dof_idx == NumLib::MeshComponentMap::nop)
94  {
95  // TODO extend Python BC to mixed FEM ansatz functions
96  OGS_FATAL(
97  "No d.o.f. found for (node={:d}, var={:d}, "
98  "comp={:d}). "
99  "That might be due to the use of mixed FEM ansatz "
100  "functions, which is currently not supported by "
101  "the implementation of Python BCs. That excludes, "
102  "e.g., the HM process.",
103  bulk_node_id, var, comp);
104  }
105  primary_variables_mat(element_node_id, global_component) =
106  (*x[process_id])[dof_idx];
107  }
108  }
109  }
110 
111  Eigen::VectorXd local_rhs = Eigen::VectorXd::Zero(num_nodes);
112  Eigen::MatrixXd local_Jac =
113  Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
114  std::vector<double> prim_vars_data(num_comp_total);
115  auto prim_vars = MathLib::toVector(prim_vars_data);
116 
117  for (unsigned ip = 0; ip < num_integration_points; ip++)
118  {
119  auto const& N = Base::_ns_and_weights[ip].N;
120  auto const& w = Base::_ns_and_weights[ip].weight;
121  auto const coords = fe.interpolateCoordinates(N);
122  prim_vars =
123  N * primary_variables_mat; // Assumption: all primary variables
124  // have same shape functions.
125  auto const flag_flux_dFlux =
126  _data.bc_object->getFlux(t, coords, prim_vars_data);
128  {
129  // getFlux() is not overridden in Python, so we can skip the
130  // whole BC assembly (i.e., for all boundary elements).
131  throw MethodNotOverriddenInDerivedClassException{};
132  }
133 
134  if (!std::get<0>(flag_flux_dFlux))
135  {
136  // No flux value for this integration point. Skip assembly of
137  // the entire element.
138  return;
139  }
140  auto const flux = std::get<1>(flag_flux_dFlux);
141  auto const& dFlux = std::get<2>(flag_flux_dFlux);
142 
143  local_rhs.noalias() += N * (flux * w);
144 
145  if (static_cast<int>(dFlux.size()) != num_comp_total)
146  {
147  // This strict check is technically mandatory only if a Jacobian
148  // is assembled. However, it is done as a consistency check also
149  // for cases without Jacobian assembly.
150  OGS_FATAL(
151  "The Python BC must return the derivative of the flux "
152  "w.r.t. each primary variable. {:d} components expected. "
153  "{:d} components returned from Python.",
154  num_comp_total, dFlux.size());
155  }
156 
157  if (Jac)
158  {
159  for (int comp = 0; comp < num_comp_total; ++comp)
160  {
161  auto const top = 0;
162  auto const left = comp * num_nodes;
163  auto const width = num_nodes;
164  auto const height = num_nodes;
165  // The assignment -= takes into account the sign convention
166  // of 1st-order in time ODE systems in OpenGeoSys.
167  local_Jac.block(top, left, width, height).noalias() -=
168  N.transpose() * (dFlux[comp] * w) * N;
169  }
170  }
171  }
172 
173  auto const& indices_specific_component =
174  dof_table_boundary(boundary_element_id, _data.global_component_id)
175  .rows;
176  b.add(indices_specific_component, local_rhs);
177 
178  if (Jac)
179  {
180  // only assemble a block of the Jacobian, not the whole local matrix
181  auto const indices_all_components =
182  NumLib::getIndices(boundary_element_id, dof_table_boundary);
184  indices_specific_component, indices_all_components};
185  Jac->add(rci, local_Jac);
186  }
187  }
#define OGS_FATAL(...)
Definition: Error.h:26
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:87
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:80
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual const Node * getNode(unsigned idx) const =0
virtual unsigned getNumberOfNodes() const =0
Properties & getProperties()
Definition: Mesh.h:123
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
int getNumberOfVariableComponents(int variable_id) const
static NUMLIB_EXPORT GlobalIndexType const nop
std::vector< NAndWeight, Eigen::aligned_allocator< NAndWeight > > const _ns_and_weights
virtual std::tuple< bool, double, std::vector< double > > getFlux(double, std::array< double, 3 >, std::vector< double > const &) const
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
NumLib::TemplateIsoparametric< ShapeFunction, ShapeMatricesType > createIsoparametricFiniteElement(MeshLib::Element const &e)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
NumLib::LocalToGlobalIndexMap const & dof_table_bulk
DOF table of the entire domain.
PythonBoundaryConditionPythonSideInterface * bc_object
Python object computing BC values.
const MeshLib::Mesh & boundary_mesh
The boundary mesh, i.e., the domain of this BC.
std::size_t const bulk_mesh_id
Mesh ID of the entire domain.

References ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::_data, ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::_element, ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::_integration_method, ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::_ns_and_weights, MathLib::EigenMatrix::add(), MathLib::EigenVector::add(), ProcessLib::PythonBoundaryConditionData::bc_object, ProcessLib::PythonBoundaryConditionData::boundary_mesh, ProcessLib::PythonBoundaryConditionData::bulk_mesh_id, NumLib::createIsoparametricFiniteElement(), ProcessLib::PythonBoundaryConditionData::dof_table_bulk, ProcessLib::PythonBoundaryConditionPythonSideInterface::getFlux(), NumLib::LocalToGlobalIndexMap::getGlobalComponent(), NumLib::LocalToGlobalIndexMap::getGlobalIndex(), MathLib::Point3dWithID::getID(), NumLib::getIndices(), MeshLib::Element::getNode(), NumLib::LocalToGlobalIndexMap::getNumberOfGlobalComponents(), MeshLib::Element::getNumberOfNodes(), NumLib::LocalToGlobalIndexMap::getNumberOfVariableComponents(), NumLib::LocalToGlobalIndexMap::getNumberOfVariables(), MeshLib::Mesh::getProperties(), ProcessLib::PythonBoundaryConditionData::global_component_id, ProcessLib::PythonBoundaryConditionPythonSideInterface::isOverriddenNatural(), MeshLib::Node, NumLib::MeshComponentMap::nop, OGS_FATAL, and MathLib::toVector().

Member Data Documentation

◆ _data

template<typename ShapeFunction , typename IntegrationMethod , int GlobalDim>
PythonBoundaryConditionData const& ProcessLib::PythonBoundaryConditionLocalAssembler< ShapeFunction, IntegrationMethod, GlobalDim >::_data
private

The documentation for this class was generated from the following file: