OGS 6.2.1-97-g73d1aeda3
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 =
96  ShapeFunction, ShapeMatricesType>(_element);
97 
98  unsigned const num_integration_points =
99  _integration_method.getNumberOfPoints();
100  auto const num_var = dof_table_source_term.getNumberOfVariables();
101  auto const num_nodes = ShapeFunction::NPOINTS;
102  auto const num_comp_total =
103  dof_table_source_term.getNumberOfComponents();
104 
105  auto const& bulk_node_ids_map =
106  *_data.source_term_mesh.getProperties()
107  .template getPropertyVector<std::size_t>("bulk_node_ids");
108 
109  // gather primary variables
110  typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
111  Eigen::Dynamic>
112  primary_variables_mat(num_nodes, num_comp_total);
113  for (int var = 0; var < num_var; ++var)
114  {
115  auto const num_comp =
116  dof_table_source_term.getNumberOfVariableComponents(var);
117  for (int comp = 0; comp < num_comp; ++comp)
118  {
119  auto const global_component =
120  dof_table_source_term.getGlobalComponent(var, comp);
121 
122  for (unsigned element_node_id = 0; element_node_id < num_nodes;
123  ++element_node_id)
124  {
125  auto const* const node = _element.getNode(element_node_id);
126  auto const boundary_node_id = node->getID();
127  auto const bulk_node_id =
128  bulk_node_ids_map[boundary_node_id];
129  MeshLib::Location loc{_data.bulk_mesh_id,
131  bulk_node_id};
132  auto const dof_idx =
133  dof_table_source_term.getGlobalIndex(loc, var, comp);
134  if (dof_idx == NumLib::MeshComponentMap::nop)
135  {
136  // TODO extend Python BC to mixed FEM ansatz functions
137  OGS_FATAL(
138  "No d.o.f. found for (node=%d, var=%d, comp=%d). "
139  "That might be due to the use of mixed FEM ansatz "
140  "functions, which is currently not supported by "
141  "the implementation of Python BCs. That excludes, "
142  "e.g., the HM process.",
143  bulk_node_id, var, comp);
144  }
145  primary_variables_mat(element_node_id, global_component) =
146  x[dof_idx];
147  }
148  }
149  }
150 
151  NodalRowVectorType local_rhs = Eigen::VectorXd::Zero(num_nodes);
152  NodalMatrixType local_Jac =
153  Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
154  std::vector<double> prim_vars_data(num_comp_total);
155  auto prim_vars = MathLib::toVector(prim_vars_data);
156 
157  for (unsigned ip = 0; ip < num_integration_points; ip++)
158  {
159  auto const& ip_data = _ip_data[ip];
160  auto const& N = ip_data.N;
161  auto const& coords = fe.interpolateCoordinates(N);
162  // Assumption: all primary variables have same shape functions.
163  prim_vars = N * primary_variables_mat;
164 
165  auto const& w = ip_data.integration_weight;
166  auto const flux_dflux =
167  _data.source_term_object->getFlux(t, coords, prim_vars_data);
168  auto const& flux = flux_dflux.first;
169  auto const& dflux = flux_dflux.second;
170  local_rhs.noalias() += N * (flux * w);
171 
172  if (static_cast<int>(dflux.size()) != num_comp_total)
173  {
174  // This strict check is technically mandatory only if a
175  // Jacobian is assembled. However, it is done as a
176  // consistency check also for cases without Jacobian
177  // assembly.
178  OGS_FATAL(
179  "The Python source term must return the derivative of "
180  "the flux w.r.t. each primary variable. %d components "
181  "expected. %d components returned from Python.",
182  num_comp_total, dflux.size());
183  }
184 
185  if (Jac)
186  {
187  for (int comp = 0; comp < num_comp_total; ++comp)
188  {
189  auto const top = 0;
190  auto const left = comp * num_nodes;
191  auto const width = num_nodes;
192  auto const height = num_nodes;
193  // The assignement -= takes into account the sign convention
194  // of 1st-order in time ODE systems in OpenGeoSys.
195  local_Jac.block(top, left, width, height).noalias() -=
196  ip_data.N.transpose() * (dflux[comp] * w) * N;
197  }
198  }
199  }
200 
201  auto const& indices_specific_component =
202  dof_table_source_term(source_term_element_id,
203  _data.global_component_id)
204  .rows;
205  b.add(indices_specific_component, local_rhs);
206 
207  if (Jac)
208  {
209  // only assemble a block of the Jacobian, not the whole local matrix
210  auto const indices_all_components = NumLib::getIndices(
211  source_term_element_id, dof_table_source_term);
213  indices_specific_component, indices_all_components};
214  Jac->add(rci, local_Jac);
215  }
216  }
217 
218 private:
222 
223  IntegrationMethod const _integration_method;
224  std::vector<
226  Eigen::aligned_allocator<IntegrationPointData<NodalRowVectorType>>>
228 };
229 
230 } // namespace Python
231 } // namespace SourceTerms
232 } // namespace ProcessLib
int getGlobalComponent(int const variable_id, int const component_id) const
NumLib::TemplateIsoparametric< ShapeFunction, ShapeMatricesType > createIsoparametricFiniteElement(MeshLib::Element const &e)
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
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