OGS 6.2.0-97-g4a610c866
PythonSourceTermLocalAssembler.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include "PythonSourceTerm.h"
13 
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,
50  unsigned GlobalDim>
53 {
54 public:
58 
60  MeshLib::Element const& e,
61  std::size_t const /*local_matrix_size*/,
62  bool is_axially_symmetric,
63  unsigned const integration_order,
64  PythonSourceTermData const& data)
65  : _data(data),
66  _element(e),
67  _is_axially_symmetric(is_axially_symmetric),
68  _integration_method(integration_order)
69  {
70  unsigned const n_integration_points =
71  _integration_method.getNumberOfPoints();
72 
73  auto const shape_matrices =
74  initShapeMatrices<ShapeFunction, ShapeMatricesType,
75  IntegrationMethod, GlobalDim>(
76  _element, _is_axially_symmetric, _integration_method);
77 
78  for (unsigned ip = 0; ip < n_integration_points; ip++)
79  {
80  _ip_data.emplace_back(
81  shape_matrices[ip].N,
82  _integration_method.getWeightedPoint(ip).getWeight() *
83  shape_matrices[ip].integralMeasure *
84  shape_matrices[ip].detJ);
85  }
86  }
87 
88  void assemble(std::size_t const source_term_element_id,
89  NumLib::LocalToGlobalIndexMap const& dof_table_source_term,
90  double const t, const GlobalVector& x, GlobalVector& b,
91  GlobalMatrix* Jac) override
92  {
93  using ShapeMatricesType =
95  using FemType =
97  FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
98  &_element));
99 
100  unsigned const num_integration_points =
101  _integration_method.getNumberOfPoints();
102  auto const num_var = dof_table_source_term.getNumberOfVariables();
103  auto const num_nodes = ShapeFunction::NPOINTS;
104  auto const num_comp_total =
105  dof_table_source_term.getNumberOfComponents();
106 
107  auto const& bulk_node_ids_map =
108  *_data.source_term_mesh.getProperties()
109  .template getPropertyVector<std::size_t>("bulk_node_ids");
110 
111  // gather primary variables
112  typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
113  Eigen::Dynamic>
114  primary_variables_mat(num_nodes, num_comp_total);
115  for (int var = 0; var < num_var; ++var)
116  {
117  auto const num_comp =
118  dof_table_source_term.getNumberOfVariableComponents(var);
119  for (int comp = 0; comp < num_comp; ++comp)
120  {
121  auto const global_component =
122  dof_table_source_term.getGlobalComponent(var, comp);
123 
124  for (unsigned element_node_id = 0; element_node_id < num_nodes;
125  ++element_node_id)
126  {
127  auto const* const node = _element.getNode(element_node_id);
128  auto const boundary_node_id = node->getID();
129  auto const bulk_node_id =
130  bulk_node_ids_map[boundary_node_id];
131  MeshLib::Location loc{_data.bulk_mesh_id,
133  bulk_node_id};
134  auto const dof_idx =
135  dof_table_source_term.getGlobalIndex(loc, var, comp);
136  if (dof_idx == NumLib::MeshComponentMap::nop)
137  {
138  // TODO extend Python BC to mixed FEM ansatz functions
139  OGS_FATAL(
140  "No d.o.f. found for (node=%d, var=%d, comp=%d). "
141  "That might be due to the use of mixed FEM ansatz "
142  "functions, which is currently not supported by "
143  "the implementation of Python BCs. That excludes, "
144  "e.g., the HM process.",
145  bulk_node_id, var, comp);
146  }
147  primary_variables_mat(element_node_id, global_component) =
148  x[dof_idx];
149  }
150  }
151  }
152 
153  NodalRowVectorType local_rhs = Eigen::VectorXd::Zero(num_nodes);
154  NodalMatrixType local_Jac =
155  Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
156  std::vector<double> prim_vars_data(num_comp_total);
157  auto prim_vars = MathLib::toVector(prim_vars_data);
158 
159  for (unsigned ip = 0; ip < num_integration_points; ip++)
160  {
161  auto const& ip_data = _ip_data[ip];
162  auto const& N = ip_data.N;
163  auto const& coords = fe.interpolateCoordinates(N);
164  // Assumption: all primary variables have same shape functions.
165  prim_vars = N * primary_variables_mat;
166 
167  auto const& w = ip_data.integration_weight;
168  auto const flux_dflux =
169  _data.source_term_object->getFlux(t, coords, prim_vars_data);
170  auto const& flux = flux_dflux.first;
171  auto const& dflux = flux_dflux.second;
172  local_rhs.noalias() += N * (flux * w);
173 
174  if (static_cast<int>(dflux.size()) != num_comp_total)
175  {
176  // This strict check is technically mandatory only if a
177  // Jacobian is assembled. However, it is done as a
178  // consistency check also for cases without Jacobian
179  // assembly.
180  OGS_FATAL(
181  "The Python source term must return the derivative of "
182  "the flux w.r.t. each primary variable. %d components "
183  "expected. %d components returned from Python.",
184  num_comp_total, dflux.size());
185  }
186 
187  if (Jac)
188  {
189  for (int comp = 0; comp < num_comp_total; ++comp)
190  {
191  auto const top = 0;
192  auto const left = comp * num_nodes;
193  auto const width = num_nodes;
194  auto const height = num_nodes;
195  // The assignement -= takes into account the sign convention
196  // of 1st-order in time ODE systems in OpenGeoSys.
197  local_Jac.block(top, left, width, height).noalias() -=
198  ip_data.N.transpose() * (dflux[comp] * w) * N;
199  }
200  }
201  }
202 
203  auto const& indices_specific_component =
204  dof_table_source_term(source_term_element_id,
205  _data.global_component_id)
206  .rows;
207  b.add(indices_specific_component, local_rhs);
208 
209  if (Jac)
210  {
211  // only assemble a block of the Jacobian, not the whole local matrix
212  auto const indices_all_components = NumLib::getIndices(
213  source_term_element_id, dof_table_source_term);
215  indices_specific_component, indices_all_components};
216  Jac->add(rci, local_Jac);
217  }
218  }
219 
220 private:
224 
225  IntegrationMethod const _integration_method;
226  std::vector<
228  Eigen::aligned_allocator<IntegrationPointData<NodalRowVectorType>>>
230 };
231 
232 } // namespace Python
233 } // namespace SourceTerms
234 } // namespace ProcessLib
int getGlobalComponent(int const variable_id, int const component_id) const
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
int getNumberOfVariableComponents(int variable_id) const
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< IntegrationPointData< NodalRowVectorType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType > > > _ip_data
Groups data used by source terms, in particular by the local assemblers.
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
PythonSourceTermLocalAssembler(MeshLib::Element const &e, std::size_t const, bool is_axially_symmetric, unsigned const integration_order, PythonSourceTermData const &data)
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Template class for isoparametric elements.
IntegrationPointData(NodalRowVectorType const &N_, double const &integration_weight_)
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
static NUMLIB_EXPORT GlobalIndexType const nop