OGS 6.1.0-1721-g6382411ad
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 // /// Lines ///////////////////////////////////
175 #if OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
176  _builder[std::type_index(typeid(MeshLib::Line3))] =
177  makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
178 #endif
179  }
180 
186  void operator()(std::size_t const id,
187  MeshLib::Element const& mesh_item,
188  LADataIntfPtr& data_ptr,
189  ConstructorArgs&&... args) const
190  {
191  auto const type_idx = std::type_index(typeid(mesh_item));
192  auto const it = _builder.find(type_idx);
193 
194  if (it == _builder.end())
195  {
196  OGS_FATAL(
197  "You are trying to build a local assembler for an unknown mesh "
198  "element type (%s)."
199  " Maybe you have disabled this mesh element type in your build "
200  "configuration or this process requires higher order elements.",
201  type_idx.name());
202  }
203 
204  auto const n_local_dof = _dof_table.getNumberOfElementDOF(id);
205  auto const varIDs = _dof_table.getElementVariableIDs(id);
206  bool const isPressureDeactivated = (varIDs.front() != 0);
207  std::vector<int> involved_varIDs; // including deactived elements
208  involved_varIDs.reserve(varIDs.size() + 1);
209  if (isPressureDeactivated)
210  {
211  involved_varIDs.push_back(0); // always pressure come in
212  }
213  involved_varIDs.insert(involved_varIDs.end(), varIDs.begin(),
214  varIDs.end());
215 
216  std::vector<unsigned> dofIndex_to_localIndex;
217 
218  // matrix and fracture assemblers with enrichments
219  dofIndex_to_localIndex.resize(n_local_dof);
220  std::vector<unsigned> vec_n_element_nodes;
221  // TODO how to get the shape function order for each variable?
222  vec_n_element_nodes.push_back(
223  mesh_item.getNumberOfBaseNodes()); // pressure
224  auto const max_varID = *std::max_element(varIDs.begin(), varIDs.end());
225  for (int i = 1; i < max_varID + 1; i++)
226  {
227  vec_n_element_nodes.push_back(
228  mesh_item.getNumberOfNodes()); // displacements
229  }
230 
231  unsigned local_id = 0;
232  unsigned dof_id = 0;
233  for (unsigned i = 0; i < involved_varIDs.size(); i++)
234  {
235  auto const var_id = involved_varIDs[i];
236  auto const n_var_comp =
238  auto const n_var_element_nodes = vec_n_element_nodes[i];
239  for (int var_comp_id = 0; var_comp_id < n_var_comp; var_comp_id++)
240  {
241  auto const& ms = _dof_table.getMeshSubset(var_id, var_comp_id);
242  auto const mesh_id = ms.getMeshID();
243  for (unsigned k = 0; k < n_var_element_nodes; k++)
244  {
245  MeshLib::Location l(mesh_id,
247  mesh_item.getNodeIndex(k));
248  auto global_index =
249  _dof_table.getGlobalIndex(l, var_id, var_comp_id);
250  if (global_index != NumLib::MeshComponentMap::nop)
251  {
252  dofIndex_to_localIndex[dof_id++] = local_id;
253  }
254  local_id++;
255  }
256  }
257  }
258 
259  data_ptr = it->second(mesh_item, involved_varIDs.size(), n_local_dof,
260  dofIndex_to_localIndex,
261  std::forward<ConstructorArgs>(args)...);
262  }
263 
264 private:
265  using LADataBuilder = std::function<LADataIntfPtr(
266  MeshLib::Element const& e,
267  std::size_t const n_variables,
268  std::size_t const local_matrix_size,
269  std::vector<unsigned> const& dofIndex_to_localIndex,
270  ConstructorArgs&&...)>;
271 
272  template <typename ShapeFunctionDisplacement>
274  typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
275 
276  template <typename ShapeFunctionDisplacement,
277  typename ShapeFunctionPressure>
278  using LADataMatrix = LocalAssemblerDataMatrix<
279  ShapeFunctionDisplacement, ShapeFunctionPressure,
281 
282  template <typename ShapeFunctionDisplacement,
283  typename ShapeFunctionPressure>
284  using LADataMatrixNearFracture = LocalAssemblerDataMatrixNearFracture<
285  ShapeFunctionDisplacement, ShapeFunctionPressure,
287 
288  template <typename ShapeFunctionDisplacement,
289  typename ShapeFunctionPressure>
290  using LAFractureData = LocalAssemblerDataFracture<
291  ShapeFunctionDisplacement, ShapeFunctionPressure,
293 
297  template <typename ShapeFunctionDisplacement>
299  {
300  return makeLocalAssemblerBuilder<ShapeFunctionDisplacement>(
301  static_cast<std::integral_constant<
302  bool, (GlobalDim >= ShapeFunctionDisplacement::DIM)>*>(
303  nullptr));
304  }
305 
307  std::unordered_map<std::type_index, LADataBuilder> _builder;
308 
310 
311 private:
318  template <typename ShapeFunctionDisplacement>
319  static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
320  {
321  // (Lower order elements = Order(ShapeFunctionDisplacement) - 1).
322  using ShapeFunctionPressure =
324  return [](MeshLib::Element const& e,
325  std::size_t const n_variables,
326  std::size_t const local_matrix_size,
327  std::vector<unsigned> const& dofIndex_to_localIndex,
328  ConstructorArgs&&... args) {
329  if (e.getDimension() == GlobalDim)
330  {
331  if (n_variables == 2)
332  {
333  return LADataIntfPtr{
334  new LADataMatrix<ShapeFunctionDisplacement,
335  ShapeFunctionPressure>{
336  e, n_variables, local_matrix_size,
337  dofIndex_to_localIndex,
338  std::forward<ConstructorArgs>(args)...}};
339  }
340  return LADataIntfPtr{new LADataMatrixNearFracture<
341  ShapeFunctionDisplacement, ShapeFunctionPressure>{
342  e, n_variables, local_matrix_size, dofIndex_to_localIndex,
343  std::forward<ConstructorArgs>(args)...}};
344  }
345  return LADataIntfPtr{new LAFractureData<ShapeFunctionDisplacement,
346  ShapeFunctionPressure>{
347  e, local_matrix_size, dofIndex_to_localIndex,
348  std::forward<ConstructorArgs>(args)...}};
349  };
350  }
351 
354  template <typename ShapeFunctionDisplacement>
355  static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
356  {
357  return nullptr;
358  }
359 };
360 
361 } // namespace HydroMechanics
362 } // namespace LIE
363 } // namespace ProcessLib
364 
365 #undef ENABLED_ELEMENT_TYPE_SIMPLEX
366 #undef ENABLED_ELEMENT_TYPE_CUBOID
367 #undef ENABLED_ELEMENT_TYPE_PYRAMID
368 #undef ENABLED_ELEMENT_TYPE_PRISM
369 #undef ENABLED_ELEMENT_TYPE_TRI
370 #undef ENABLED_ELEMENT_TYPE_QUAD
371 #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:164
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:71
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 *)