OGS
PythonSourceTermLocalAssembler.h
Go to the documentation of this file.
1 
11 #pragma once
12 
18 #include "PythonSourceTerm.h"
21 
22 namespace ProcessLib
23 {
24 namespace SourceTerms
25 {
26 namespace Python
27 {
31 {
32 };
33 
34 template <typename NodalRowVectorType>
36 {
37  IntegrationPointData(NodalRowVectorType const& N_,
38  double const& integration_weight_)
39  : N(N_), integration_weight(integration_weight_)
40  {
41  }
42 
43  NodalRowVectorType const N;
44  double const integration_weight;
45 
47 };
48 
49 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
52 {
53 public:
57 
59  MeshLib::Element const& e,
60  std::size_t const /*local_matrix_size*/,
61  bool is_axially_symmetric,
62  unsigned const integration_order,
63  PythonSourceTermData const& data)
64  : _data(data),
65  _element(e),
66  _is_axially_symmetric(is_axially_symmetric),
67  _integration_method(integration_order)
68  {
69  unsigned const n_integration_points =
70  _integration_method.getNumberOfPoints();
71 
72  auto const shape_matrices =
74  GlobalDim>(
76 
77  for (unsigned ip = 0; ip < n_integration_points; ip++)
78  {
79  _ip_data.emplace_back(
80  shape_matrices[ip].N,
81  _integration_method.getWeightedPoint(ip).getWeight() *
82  shape_matrices[ip].integralMeasure *
83  shape_matrices[ip].detJ);
84  }
85  }
86 
87  void assemble(std::size_t const source_term_element_id,
88  NumLib::LocalToGlobalIndexMap const& dof_table_source_term,
89  double const t, const GlobalVector& x, GlobalVector& b,
90  GlobalMatrix* Jac) override
91  {
92  using ShapeMatricesType =
95  ShapeFunction, ShapeMatricesType>(_element);
96 
97  unsigned const num_integration_points =
98  _integration_method.getNumberOfPoints();
99  auto const num_var = dof_table_source_term.getNumberOfVariables();
100  auto const num_nodes = ShapeFunction::NPOINTS;
101  auto const num_comp_total =
102  dof_table_source_term.getNumberOfGlobalComponents();
103 
104  // gather primary variables
105  typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
106  Eigen::Dynamic>
107  primary_variables_mat(num_nodes, num_comp_total);
108  for (int var = 0; var < num_var; ++var)
109  {
110  auto const num_comp =
111  dof_table_source_term.getNumberOfVariableComponents(var);
112  for (int comp = 0; comp < num_comp; ++comp)
113  {
114  auto const global_component =
115  dof_table_source_term.getGlobalComponent(var, comp);
116 
117  for (unsigned element_node_id = 0; element_node_id < num_nodes;
118  ++element_node_id)
119  {
120  auto const boundary_node_id =
121  _element.getNode(element_node_id)->getID();
124  boundary_node_id};
125  auto const dof_idx =
126  dof_table_source_term.getGlobalIndex(loc, var, comp);
127  if (dof_idx == NumLib::MeshComponentMap::nop)
128  {
129  // TODO extend Python BC to mixed FEM ansatz functions
130  OGS_FATAL(
131  "No d.o.f. found for (node={:d}, var={:d}, "
132  "comp={:d}). "
133  "That might be due to the use of mixed FEM ansatz "
134  "functions, which is currently not supported by "
135  "the implementation of Python BCs. That excludes, "
136  "e.g., the HM process.",
137  boundary_node_id, var, comp);
138  }
139  primary_variables_mat(element_node_id, global_component) =
140  x[dof_idx];
141  }
142  }
143  }
144 
145  NodalRowVectorType local_rhs = Eigen::VectorXd::Zero(num_nodes);
146  NodalMatrixType local_Jac =
147  Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
148  std::vector<double> prim_vars_data(num_comp_total);
149  auto prim_vars = MathLib::toVector(prim_vars_data);
150 
151  for (unsigned ip = 0; ip < num_integration_points; ip++)
152  {
153  auto const& ip_data = _ip_data[ip];
154  auto const& N = ip_data.N;
155  auto const& coords = fe.interpolateCoordinates(N);
156  // Assumption: all primary variables have same shape functions.
157  prim_vars = N * primary_variables_mat;
158 
159  auto const& w = ip_data.integration_weight;
160  auto const flux_dflux =
161  _data.source_term_object->getFlux(t, coords, prim_vars_data);
162  auto const& flux = flux_dflux.first;
163  auto const& dflux = flux_dflux.second;
164  local_rhs.noalias() += N * (flux * w);
165 
166  if (static_cast<int>(dflux.size()) != num_comp_total)
167  {
168  // This strict check is technically mandatory only if a
169  // Jacobian is assembled. However, it is done as a
170  // consistency check also for cases without Jacobian
171  // assembly.
172  OGS_FATAL(
173  "The Python source term must return the derivative of "
174  "the flux w.r.t. each primary variable. {:d} components "
175  "expected. {:d} components returned from Python.",
176  num_comp_total, dflux.size());
177  }
178 
179  if (Jac)
180  {
181  for (int comp = 0; comp < num_comp_total; ++comp)
182  {
183  auto const top = 0;
184  auto const left = comp * num_nodes;
185  auto const width = num_nodes;
186  auto const height = num_nodes;
187  // The assignment -= takes into account the sign convention
188  // of 1st-order in time ODE systems in OpenGeoSys.
189  local_Jac.block(top, left, width, height).noalias() -=
190  ip_data.N.transpose() * (dflux[comp] * w) * N;
191  }
192  }
193  }
194 
195  auto const& indices_specific_component =
196  dof_table_source_term(source_term_element_id,
198  .rows;
199  b.add(indices_specific_component, local_rhs);
200 
201  if (Jac)
202  {
203  // only assemble a block of the Jacobian, not the whole local matrix
204  auto const indices_all_components = NumLib::getIndices(
205  source_term_element_id, dof_table_source_term);
207  indices_specific_component, indices_all_components};
208  Jac->add(rci, local_Jac);
209  }
210  }
211 
212 private:
216 
217  IntegrationMethod const _integration_method;
218  std::vector<
220  Eigen::aligned_allocator<IntegrationPointData<NodalRowVectorType>>>
222 };
223 
224 } // namespace Python
225 } // namespace SourceTerms
226 } // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:87
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:79
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual const Node * getNode(unsigned idx) const =0
int getGlobalComponent(int const variable_id, int const component_id) const
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< IntegrationPointData< NodalRowVectorType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType > > > _ip_data
void assemble(std::size_t const source_term_element_id, NumLib::LocalToGlobalIndexMap const &dof_table_source_term, double const t, const GlobalVector &x, GlobalVector &b, GlobalMatrix *Jac) override
PythonSourceTermLocalAssembler(MeshLib::Element const &e, std::size_t const, bool is_axially_symmetric, unsigned const integration_order, PythonSourceTermData const &data)
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
virtual std::pair< double, std::vector< double > > getFlux(double, std::array< double, 3 > const &, std::vector< double > const &) const =0
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
static const double t
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)
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)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
IntegrationPointData(NodalRowVectorType const &N_, double const &integration_weight_)
Groups data used by source terms, in particular by the local assemblers.
std::size_t const source_term_mesh_id
Mesh ID of the entire domain.
PythonSourceTermPythonSideInterface * source_term_object
Python object computing source term values.