OGS 6.2.1-97-g73d1aeda3
LocalDataInitializer.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <functional>
13 #include <memory>
14 #include <type_traits>
15 #include <typeindex>
16 #include <typeinfo>
17 #include <unordered_map>
18 
23 
24 #ifndef OGS_MAX_ELEMENT_DIM
25 static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
26 #endif
27 
28 #ifndef OGS_MAX_ELEMENT_ORDER
29 static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
30 #endif
31 
32 // The following macros decide which element types will be compiled, i.e.
33 // which element types will be available for use in simulations.
34 
35 #ifdef OGS_ENABLE_ELEMENT_SIMPLEX
36 #define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
37 #else
38 #define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
39 #endif
40 
41 #ifdef OGS_ENABLE_ELEMENT_CUBOID
42 #define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
43 #else
44 #define ENABLED_ELEMENT_TYPE_CUBOID 0u
45 #endif
46 
47 #ifdef OGS_ENABLE_ELEMENT_PRISM
48 #define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
49 #else
50 #define ENABLED_ELEMENT_TYPE_PRISM 0u
51 #endif
52 
53 #ifdef OGS_ENABLE_ELEMENT_PYRAMID
54 #define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
55 #else
56 #define ENABLED_ELEMENT_TYPE_PYRAMID 0u
57 #endif
58 
59 // Dependent element types.
60 // Faces of tets, pyramids and prisms are triangles
61 #define ENABLED_ELEMENT_TYPE_TRI \
62  ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
63  (ENABLED_ELEMENT_TYPE_PRISM))
64 // Faces of hexes, pyramids and prisms are quads
65 #define ENABLED_ELEMENT_TYPE_QUAD \
66  ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
67  (ENABLED_ELEMENT_TYPE_PRISM))
68 
69 // All enabled element types
70 #define OGS_ENABLED_ELEMENTS \
71  ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
72  (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
73 
74 // Include only what is needed (Well, the conditions are not sharp).
75 #if OGS_ENABLED_ELEMENTS != 0
78 #endif
79 
80 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
83 #endif
84 
85 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
88 #endif
89 
90 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
93 #endif
94 
95 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
99 #endif
100 
101 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
104 #endif
105 
106 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
109 #endif
110 
111 namespace ProcessLib
112 {
113 namespace LIE
114 {
115 namespace HydroMechanics
116 {
126 template <typename LocalAssemblerInterface,
127  template <typename, typename, typename, int>
128  class LocalAssemblerDataMatrix,
129  template <typename, typename, typename, int>
130  class LocalAssemblerDataMatrixNearFracture,
131  template <typename, typename, typename, int>
132  class LocalAssemblerDataFracture,
133  int GlobalDim, typename... ConstructorArgs>
135 {
136 public:
137  using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
138 
140  const unsigned shapefunction_order)
141  : _dof_table(dof_table)
142  {
143  if (shapefunction_order != 2)
144  {
145  OGS_FATAL(
146  "The given shape function order %d is not supported.\nOnly "
147  "shape functions of order 2 are supported.",
148  shapefunction_order);
149  }
150  // /// Quads and Hexahedra ///////////////////////////////////
151 
152 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
153  OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
154  _builder[std::type_index(typeid(MeshLib::Quad8))] =
155  makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
156  _builder[std::type_index(typeid(MeshLib::Quad9))] =
157  makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
158 #endif
159 
160 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
161  OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
162  _builder[std::type_index(typeid(MeshLib::Hex20))] =
163  makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
164 #endif
165 
166  // /// Simplices ////////////////////////////////////////////////
167 
168 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
169  OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
170  _builder[std::type_index(typeid(MeshLib::Tri6))] =
171  makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
172 #endif
173 
174 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
175  OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
176  _builder[std::type_index(typeid(MeshLib::Tet10))] =
177  makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
178 #endif
179 
180  // /// Prisms ////////////////////////////////////////////////////
181 
182 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
183  OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
184  _builder[std::type_index(typeid(MeshLib::Prism15))] =
185  makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
186 #endif
187 
188  // /// Pyramids //////////////////////////////////////////////////
189 
190 #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
191  OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
192  _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
193  makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
194 #endif
195  // /// Lines ///////////////////////////////////
196 
197 // /// Lines ///////////////////////////////////
198 #if OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
199  _builder[std::type_index(typeid(MeshLib::Line3))] =
200  makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
201 #endif
202  }
203 
209  void operator()(std::size_t const id,
210  MeshLib::Element const& mesh_item,
211  LADataIntfPtr& data_ptr,
212  ConstructorArgs&&... args) const
213  {
214  auto const type_idx = std::type_index(typeid(mesh_item));
215  auto const it = _builder.find(type_idx);
216 
217  if (it == _builder.end())
218  {
219  OGS_FATAL(
220  "You are trying to build a local assembler for an unknown mesh "
221  "element type (%s)."
222  " Maybe you have disabled this mesh element type in your build "
223  "configuration or this process requires higher order elements.",
224  type_idx.name());
225  }
226 
227  auto const n_local_dof = _dof_table.getNumberOfElementDOF(id);
228  auto const varIDs = _dof_table.getElementVariableIDs(id);
229  bool const isPressureDeactivated = (varIDs.front() != 0);
230  std::vector<int> involved_varIDs; // including deactived elements
231  involved_varIDs.reserve(varIDs.size() + 1);
232  if (isPressureDeactivated)
233  {
234  involved_varIDs.push_back(0); // always pressure come in
235  }
236  involved_varIDs.insert(involved_varIDs.end(), varIDs.begin(),
237  varIDs.end());
238 
239  std::vector<unsigned> dofIndex_to_localIndex;
240 
241  // matrix and fracture assemblers with enrichments
242  dofIndex_to_localIndex.resize(n_local_dof);
243  std::vector<unsigned> vec_n_element_nodes;
244  // TODO how to get the shape function order for each variable?
245  vec_n_element_nodes.push_back(
246  mesh_item.getNumberOfBaseNodes()); // pressure
247  auto const max_varID = *std::max_element(varIDs.begin(), varIDs.end());
248  for (int i = 1; i < max_varID + 1; i++)
249  {
250  vec_n_element_nodes.push_back(
251  mesh_item.getNumberOfNodes()); // displacements
252  }
253 
254  unsigned local_id = 0;
255  unsigned dof_id = 0;
256  for (unsigned i = 0; i < involved_varIDs.size(); i++)
257  {
258  auto const var_id = involved_varIDs[i];
259  auto const n_var_comp =
261  auto const n_var_element_nodes = vec_n_element_nodes[i];
262  for (int var_comp_id = 0; var_comp_id < n_var_comp; var_comp_id++)
263  {
264  auto const& ms = _dof_table.getMeshSubset(var_id, var_comp_id);
265  auto const mesh_id = ms.getMeshID();
266  for (unsigned k = 0; k < n_var_element_nodes; k++)
267  {
268  MeshLib::Location l(mesh_id,
270  mesh_item.getNodeIndex(k));
271  auto global_index =
272  _dof_table.getGlobalIndex(l, var_id, var_comp_id);
273  if (global_index != NumLib::MeshComponentMap::nop)
274  {
275  dofIndex_to_localIndex[dof_id++] = local_id;
276  }
277  local_id++;
278  }
279  }
280  }
281 
282  data_ptr = it->second(mesh_item, involved_varIDs.size(), n_local_dof,
283  dofIndex_to_localIndex,
284  std::forward<ConstructorArgs>(args)...);
285  }
286 
287 private:
288  using LADataBuilder = std::function<LADataIntfPtr(
289  MeshLib::Element const& e,
290  std::size_t const n_variables,
291  std::size_t const local_matrix_size,
292  std::vector<unsigned> const& dofIndex_to_localIndex,
293  ConstructorArgs&&...)>;
294 
295  template <typename ShapeFunctionDisplacement>
297  typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
298 
299  template <typename ShapeFunctionDisplacement,
300  typename ShapeFunctionPressure>
301  using LADataMatrix = LocalAssemblerDataMatrix<
302  ShapeFunctionDisplacement, ShapeFunctionPressure,
304 
305  template <typename ShapeFunctionDisplacement,
306  typename ShapeFunctionPressure>
307  using LADataMatrixNearFracture = LocalAssemblerDataMatrixNearFracture<
308  ShapeFunctionDisplacement, ShapeFunctionPressure,
310 
311  template <typename ShapeFunctionDisplacement,
312  typename ShapeFunctionPressure>
313  using LAFractureData = LocalAssemblerDataFracture<
314  ShapeFunctionDisplacement, ShapeFunctionPressure,
316 
320  template <typename ShapeFunctionDisplacement>
322  {
323  return makeLocalAssemblerBuilder<ShapeFunctionDisplacement>(
324  static_cast<std::integral_constant<
325  bool, (GlobalDim >= ShapeFunctionDisplacement::DIM)>*>(
326  nullptr));
327  }
328 
330  std::unordered_map<std::type_index, LADataBuilder> _builder;
331 
333 
334 private:
341  template <typename ShapeFunctionDisplacement>
342  static LADataBuilder makeLocalAssemblerBuilder(std::true_type* /*unused*/)
343  {
344  // (Lower order elements = Order(ShapeFunctionDisplacement) - 1).
345  using ShapeFunctionPressure =
347  return [](MeshLib::Element const& e,
348  std::size_t const n_variables,
349  std::size_t const local_matrix_size,
350  std::vector<unsigned> const& dofIndex_to_localIndex,
351  ConstructorArgs&&... args) {
352  if (e.getDimension() == GlobalDim)
353  {
354  if (n_variables == 2)
355  {
356  return LADataIntfPtr{
357  new LADataMatrix<ShapeFunctionDisplacement,
358  ShapeFunctionPressure>{
359  e, n_variables, local_matrix_size,
360  dofIndex_to_localIndex,
361  std::forward<ConstructorArgs>(args)...}};
362  }
363  return LADataIntfPtr{new LADataMatrixNearFracture<
364  ShapeFunctionDisplacement, ShapeFunctionPressure>{
365  e, n_variables, local_matrix_size, dofIndex_to_localIndex,
366  std::forward<ConstructorArgs>(args)...}};
367  }
368  return LADataIntfPtr{new LAFractureData<ShapeFunctionDisplacement,
369  ShapeFunctionPressure>{
370  e, local_matrix_size, dofIndex_to_localIndex,
371  std::forward<ConstructorArgs>(args)...}};
372  };
373  }
374 
377  template <typename ShapeFunctionDisplacement>
378  static LADataBuilder makeLocalAssemblerBuilder(std::false_type* /*unused*/)
379  {
380  return nullptr;
381  }
382 };
383 
384 } // namespace HydroMechanics
385 } // namespace LIE
386 } // namespace ProcessLib
387 
388 #undef ENABLED_ELEMENT_TYPE_SIMPLEX
389 #undef ENABLED_ELEMENT_TYPE_CUBOID
390 #undef ENABLED_ELEMENT_TYPE_PYRAMID
391 #undef ENABLED_ELEMENT_TYPE_PRISM
392 #undef ENABLED_ELEMENT_TYPE_TRI
393 #undef ENABLED_ELEMENT_TYPE_QUAD
394 #undef OGS_ENABLED_ELEMENTS
virtual unsigned getNumberOfBaseNodes() const =0
std::vector< int > getElementVariableIDs(std::size_t const mesh_item_id) const
std::function< LADataIntfPtr(MeshLib::Element const &e, std::size_t const n_variables, std::size_t const local_matrix_size, std::vector< unsigned > const &dofIndex_to_localIndex, ConstructorArgs &&...)> LADataBuilder
std::size_t getMeshID() const
return this mesh ID
Definition: MeshSubset.h:71
std::size_t getNumberOfElementDOF(std::size_t const mesh_item_id) const
static LADataBuilder makeLocalAssemblerBuilder(std::true_type *)
LocalAssemblerDataFracture< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod< ShapeFunctionDisplacement >, GlobalDim > LAFractureData
int getNumberOfVariableComponents(int variable_id) const
void operator()(std::size_t const id, MeshLib::Element const &mesh_item, LADataIntfPtr &data_ptr, ConstructorArgs &&... args) const
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
std::unique_ptr< LocalAssemblerInterface > LADataIntfPtr
std::unordered_map< std::type_index, LADataBuilder > _builder
Mapping of element types to local assembler constructors.
LocalDataInitializer(NumLib::LocalToGlobalIndexMap const &dof_table, const unsigned shapefunction_order)
std::size_t getNodeIndex(unsigned i) const
Definition: Element.cpp:182
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
MeshLib::MeshSubset const & getMeshSubset(int const variable_id, int const component_id) const
LocalAssemblerDataMatrixNearFracture< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod< ShapeFunctionDisplacement >, GlobalDim > LADataMatrixNearFracture
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
virtual unsigned getNumberOfNodes() const =0
Returns the number of all nodes including both linear and nonlinear nodes.
LocalAssemblerDataMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod< ShapeFunctionDisplacement >, GlobalDim > LADataMatrix
typename NumLib::GaussLegendreIntegrationPolicy< typename ShapeFunctionDisplacement::MeshElement >::IntegrationMethod IntegrationMethod
static NUMLIB_EXPORT GlobalIndexType const nop
static LADataBuilder makeLocalAssemblerBuilder(std::false_type *)