OGS
SmallDeformationProcess.cpp
Go to the documentation of this file.
1
12
18#include "MeshLib/Mesh.h"
19#include "MeshLib/Properties.h"
24
25namespace ProcessLib
26{
27namespace LIE
28{
29namespace SmallDeformation
30{
31template <int DisplacementDim>
33 std::string name,
34 MeshLib::Mesh& mesh,
35 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
36 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
37 unsigned const integration_order,
38 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
39 process_variables,
41 SecondaryVariableCollection&& secondary_variables)
42 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
43 integration_order, std::move(process_variables),
44 std::move(secondary_variables)),
45 _process_data(std::move(process_data))
46{
47 INFO("[LIE/M] looking for fracture elements in the given mesh");
48 std::vector<std::pair<std::size_t, std::vector<int>>>
49 vec_branch_nodeID_matIDs;
50 std::vector<std::pair<std::size_t, std::vector<int>>>
51 vec_junction_nodeID_matIDs;
55 _vec_fracture_nodes, vec_branch_nodeID_matIDs,
56 vec_junction_nodeID_matIDs);
57
58 if (_vec_fracture_mat_IDs.size() !=
59 _process_data.fracture_properties.size())
60 {
62 "The number of the given fracture properties ({:d}) are not "
63 "consistent with the number of fracture groups in a mesh ({:d}).",
64 _process_data.fracture_properties.size(),
66 }
67
68 // create a map from a material ID to a fracture ID
69 auto max_frac_mat_id = std::max_element(_vec_fracture_mat_IDs.begin(),
71 _process_data.map_materialID_to_fractureID.resize(*max_frac_mat_id + 1);
72 for (unsigned i = 0; i < _vec_fracture_mat_IDs.size(); i++)
73 {
74 _process_data.map_materialID_to_fractureID[_vec_fracture_mat_IDs[i]] =
75 i;
76 }
77
78 // create a table of connected fracture IDs for each element
79 _process_data.vec_ele_connected_fractureIDs.resize(
80 mesh.getNumberOfElements());
81 for (unsigned i = 0; i < _vec_fracture_matrix_elements.size(); i++)
82 {
83 for (auto e : _vec_fracture_matrix_elements[i])
84 {
85 _process_data.vec_ele_connected_fractureIDs[e->getID()].push_back(
86 i);
87 }
88 }
89
90 // set fracture property
91 for (auto& fracture_prop : _process_data.fracture_properties)
92 {
93 // based on the 1st element assuming a fracture forms a straight line
95 DisplacementDim,
96 *_vec_fracture_elements[fracture_prop.fracture_id][0],
97 fracture_prop);
98 }
99
100 // set branches
101 for (auto const& [vec_branch_nodeID, matID] : vec_branch_nodeID_matIDs)
102 {
103 auto const master_matId = matID[0];
104 auto const slave_matId = matID[1];
105 auto& master_frac =
106 _process_data.fracture_properties
107 [_process_data.map_materialID_to_fractureID[master_matId]];
108 auto& slave_frac =
109 _process_data.fracture_properties
110 [_process_data.map_materialID_to_fractureID[slave_matId]];
111
112 master_frac.branches_master.push_back(createBranchProperty(
113 *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
114
115 slave_frac.branches_slave.push_back(createBranchProperty(
116 *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
117 }
118
119 // set junctions
120 transform(cbegin(vec_junction_nodeID_matIDs),
121 cend(vec_junction_nodeID_matIDs),
122 back_inserter(_vec_junction_nodes),
123 [&](auto& vec_junction_nodeID_matID)
124 {
125 return const_cast<MeshLib::Node*>(
126 _mesh.getNode(vec_junction_nodeID_matID.first));
127 });
128
129 for (std::size_t i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
130 {
131 auto const& material_ids = vec_junction_nodeID_matIDs[i].second;
132 assert(material_ids.size() == 2);
133 std::array<int, 2> fracture_ids{
134 {_process_data.map_materialID_to_fractureID[material_ids[0]],
135 _process_data.map_materialID_to_fractureID[material_ids[1]]}};
136
137 _process_data.junction_properties.emplace_back(
138 i, *mesh.getNode(vec_junction_nodeID_matIDs[i].first),
139 fracture_ids);
140 }
141
142 // create a table of connected junction IDs for each element
143 _process_data.vec_ele_connected_junctionIDs.resize(
144 mesh.getNumberOfElements());
145 for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
146 {
147 auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
148 for (auto id :
150 {
151 _process_data.vec_ele_connected_junctionIDs[id].push_back(i);
152 }
153 }
154
155 // create a table of junction node and connected elements
157 vec_junction_nodeID_matIDs.size());
158 for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
159 {
160 auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
161 for (auto e : mesh.getElementsConnectedToNode(*node))
162 {
164 const_cast<MeshLib::Element*>(e));
165 }
166 }
167
168 //
169 // If Neumann BCs for the displacement_jump variable are required they need
170 // special treatment because of the levelset function. The implementation
171 // exists in the version 6.1.0 (e54815cc07ee89c81f953a4955b1c788595dd725)
172 // and was removed due to lack of applications.
173 //
174
175 MeshLib::PropertyVector<int> const* material_ids(
176 mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
177 _process_data.mesh_prop_materialIDs = material_ids;
178}
179
180template <int DisplacementDim>
182{
183 //------------------------------------------------------------
184 // prepare mesh subsets to define DoFs
185 //------------------------------------------------------------
186 // for extrapolation
187 _mesh_subset_all_nodes =
188 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
189 // regular u
190 _mesh_subset_matrix_nodes =
191 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
192 // u jump
193 for (unsigned i = 0; i < _vec_fracture_nodes.size(); i++)
194 {
195 _mesh_subset_fracture_nodes.push_back(
196 std::make_unique<MeshLib::MeshSubset const>(
197 _mesh, _vec_fracture_nodes[i]));
198 }
199 // enrichment for junctions
200 _mesh_subset_junction_nodes =
201 std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_junction_nodes);
202
203 // Collect the mesh subsets in a vector.
204 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
205 std::generate_n(std::back_inserter(all_mesh_subsets), DisplacementDim,
206 [&]() { return *_mesh_subset_matrix_nodes; });
207 for (auto const& ms : _mesh_subset_fracture_nodes)
208 {
209 std::generate_n(std::back_inserter(all_mesh_subsets),
210 DisplacementDim,
211 [&]() { return *ms; });
212 }
213 std::generate_n(std::back_inserter(all_mesh_subsets),
214 DisplacementDim,
215 [&]() { return *_mesh_subset_junction_nodes; });
216
217 std::vector<int> const vec_n_components(
218 1 + _vec_fracture_mat_IDs.size() + _vec_junction_nodes.size(),
219 DisplacementDim);
220
221 std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
222 vec_var_elements.push_back(&_vec_matrix_elements);
223 for (unsigned i = 0; i < _vec_fracture_matrix_elements.size(); i++)
224 {
225 vec_var_elements.push_back(&_vec_fracture_matrix_elements[i]);
226 }
227 for (unsigned i = 0; i < _vec_junction_fracture_matrix_elements.size(); i++)
228 {
229 vec_var_elements.push_back(&_vec_junction_fracture_matrix_elements[i]);
230 }
231
232 INFO("[LIE/M] creating a DoF table");
233 _local_to_global_index_map =
234 std::make_unique<NumLib::LocalToGlobalIndexMap>(
235 std::move(all_mesh_subsets),
236 vec_n_components,
237 vec_var_elements,
239
240 DBUG("[LIE/M] created {:d} DoF", _local_to_global_index_map->size());
241}
242
243template <int DisplacementDim>
245 MeshLib::Element const& e, MeshLib::Mesh& mesh)
246{
247 Eigen::Vector3d const pt(getCenterOfGravity(e).asEigenVector3d());
248 std::vector<FractureProperty*> e_fracture_props;
249 std::unordered_map<int, int> e_fracID_to_local;
250 unsigned tmpi = 0;
251 for (auto fid : _process_data.vec_ele_connected_fractureIDs[e.getID()])
252 {
253 e_fracture_props.push_back(&_process_data.fracture_properties[fid]);
254 e_fracID_to_local.insert({fid, tmpi++});
255 }
256 std::vector<JunctionProperty*> e_junction_props;
257 std::unordered_map<int, int> e_juncID_to_local;
258 tmpi = 0;
259 for (auto fid : _process_data.vec_ele_connected_junctionIDs[e.getID()])
260 {
261 e_junction_props.push_back(&_process_data.junction_properties[fid]);
262 e_juncID_to_local.insert({fid, tmpi++});
263 }
264 std::vector<double> const levelsets(uGlobalEnrichments(
265 e_fracture_props, e_junction_props, e_fracID_to_local, pt));
266
267 auto update_levelset_property = [&](unsigned const i, int const id,
268 unsigned const levelset_idx_offset,
269 unsigned const name_offset)
270 {
271 auto levelset_property = MeshLib::getOrCreateMeshProperty<double>(
272 const_cast<MeshLib::Mesh&>(mesh),
273 "levelset" + std::to_string(id + 1 + name_offset),
275 levelset_property->resize(mesh.getNumberOfElements());
276 (*levelset_property)[e.getID()] = levelsets[i + levelset_idx_offset];
277 };
278
279 for (unsigned i = 0; i < e_fracture_props.size(); i++)
280 {
281 update_levelset_property(i, e_fracture_props[i]->fracture_id, 0, 0);
282 }
283 for (unsigned i = 0; i < e_junction_props.size(); i++)
284 {
285 update_levelset_property(i, e_junction_props[i]->junction_id,
286 e_fracture_props.size(),
287 _process_data.fracture_properties.size());
288 }
289}
290
291template <int DisplacementDim>
293 NumLib::LocalToGlobalIndexMap const& dof_table,
294 MeshLib::Mesh const& mesh,
295 unsigned const integration_order)
296{
301 mesh.getElements(), dof_table, _local_assemblers,
302 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
303 _process_data);
304
305 // TODO move the two data members somewhere else.
306 // for extrapolation of secondary variables
307 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
308 *_mesh_subset_all_nodes};
309 _local_to_global_index_map_single_component =
310 std::make_unique<NumLib::LocalToGlobalIndexMap>(
311 std::move(all_mesh_subsets_single_component),
312 // by location order is needed for output
314
315 auto add_secondary_variable = [&](std::string const& name,
316 int const num_components,
317 auto get_ip_values_function)
318 {
319 _secondary_variables.addSecondaryVariable(
320 name,
321 makeExtrapolator(num_components, getExtrapolator(),
322 _local_assemblers,
323 std::move(get_ip_values_function)));
324 };
325
326 add_secondary_variable("sigma",
328 DisplacementDim>::RowsAtCompileTime,
329 &LocalAssemblerInterface::getIntPtSigma);
330
331 add_secondary_variable("epsilon",
333 DisplacementDim>::RowsAtCompileTime,
334 &LocalAssemblerInterface::getIntPtEpsilon);
335
336 add_secondary_variable("fracture_stress", DisplacementDim,
337 &LocalAssemblerInterface::getIntPtFractureStress);
338
339 add_secondary_variable("fracture_aperture", 1,
340 &LocalAssemblerInterface::getIntPtFractureAperture);
341
342 _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
343 const_cast<MeshLib::Mesh&>(mesh), "sigma_avg",
346 DisplacementDim>::RowsAtCompileTime);
347
348 for (MeshLib::Element const* e : _mesh.getElements())
349 {
350 if (e->getDimension() < DisplacementDim)
351 {
352 continue;
353 }
354
355 updateElementLevelSets(*e, const_cast<MeshLib::Mesh&>(mesh));
356 }
357
358 _process_data.element_local_jumps =
360 const_cast<MeshLib::Mesh&>(mesh), "local_jump_w_avg",
361 MeshLib::MeshItemType::Cell, DisplacementDim);
362
363 _process_data.element_fracture_stresses =
365 const_cast<MeshLib::Mesh&>(mesh), "fracture_stress_avg",
366 MeshLib::MeshItemType::Cell, DisplacementDim);
367
369 const_cast<MeshLib::Mesh&>(mesh), "fracture_aperture_avg",
371
372 mesh_prop_b->resize(mesh.getNumberOfElements());
373 auto const& mesh_prop_matid = *_process_data.mesh_prop_materialIDs;
374 for (auto const& fracture_prop : _process_data.fracture_properties)
375 {
376 for (MeshLib::Element const* e : _mesh.getElements())
377 {
378 if (e->getDimension() == DisplacementDim)
379 {
380 continue;
381 }
382 if (mesh_prop_matid[e->getID()] != fracture_prop.mat_id)
383 {
384 continue;
385 }
386 // Mean value for the element. This allows usage of node based
387 // properties for aperture.
388 (*mesh_prop_b)[e->getID()] =
389 fracture_prop.aperture0
390 .getNodalValuesOnElement(*e, /*time independent*/ 0)
391 .mean();
392 }
393 }
394 _process_data.mesh_prop_b = mesh_prop_b;
395}
396
397template <int DisplacementDim>
399 double const t, double const dt, std::vector<GlobalVector*> const& x,
400 GlobalVector const& x_prev, int const process_id)
401{
402 DBUG("Compute the secondary variables for SmallDeformationProcess.");
403 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
404 dof_tables.reserve(x.size());
405 std::generate_n(std::back_inserter(dof_tables), x.size(),
406 [&]() { return _local_to_global_index_map.get(); });
407
410 _local_assemblers, getActiveElementIDs(), dof_tables, t, dt, x, x_prev,
411 process_id);
412}
413
414template <int DisplacementDim>
416{
417 return false;
418}
419
420template <int DisplacementDim>
422 const double t, double const dt, std::vector<GlobalVector*> const& x,
423 std::vector<GlobalVector*> const& x_prev, int const process_id,
425{
426 DBUG("Assemble SmallDeformationProcess.");
427
428 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
429 _local_to_global_index_map.get()};
430 // Call global assembler for each local assembly item.
432 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
433 getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, &M, &K,
434 &b);
435}
436
437template <int DisplacementDim>
440 const double t, double const dt, std::vector<GlobalVector*> const& x,
441 std::vector<GlobalVector*> const& x_prev, int const process_id,
442 GlobalVector& b, GlobalMatrix& Jac)
443{
444 DBUG("AssembleWithJacobian SmallDeformationProcess.");
445
446 // Call global assembler for each local assembly item.
447 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
448 _local_to_global_index_map.get()};
451 _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x, x_prev,
452 process_id, &b, &Jac);
453}
454
455template <int DisplacementDim>
457 std::vector<GlobalVector*> const& x, double const t, double const dt,
458 const int process_id)
459{
460 DBUG("PreTimestep SmallDeformationProcess.");
461
464 _local_assemblers, getActiveElementIDs(), *_local_to_global_index_map,
465 *x[process_id], t, dt);
466}
467
468// ------------------------------------------------------------------------------------
469// template instantiation
470// ------------------------------------------------------------------------------------
471template class SmallDeformationProcess<2>;
472template class SmallDeformationProcess<3>;
473
474} // namespace SmallDeformation
475} // namespace LIE
476} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the class Properties that implements a container of properties.
Definition of the Mesh class.
Global vector based on Eigen vector.
Definition EigenVector.h:26
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
bool isAxiallySymmetric() const
Definition Mesh.h:139
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition Mesh.h:93
Properties & getProperties()
Definition Mesh.h:136
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:257
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
PropertyVector< T > const * getPropertyVector(std::string_view name) const
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
std::vector< std::vector< MeshLib::Element * > > _vec_junction_fracture_matrix_elements
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id) override
std::vector< std::vector< MeshLib::Node * > > _vec_fracture_nodes
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) override
std::vector< std::vector< MeshLib::Element * > > _vec_fracture_elements
SmallDeformationProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, SmallDeformationProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables)
void updateElementLevelSets(MeshLib::Element const &e, MeshLib::Mesh &mesh)
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
std::vector< std::vector< MeshLib::Element * > > _vec_fracture_matrix_elements
SmallDeformationProcessData< DisplacementDim > _process_data
virtual void computeSecondaryVariable(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
MeshLib::Mesh & _mesh
Definition Process.h:365
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
@ BY_LOCATION
Ordering data by spatial location.
@ BY_COMPONENT
Ordering data by component type.
void createLocalAssemblers(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
void setFractureProperty(int const dim, MeshLib::Element const &e, FractureProperty &frac_prop)
void getFractureMatrixDataInMesh(MeshLib::Mesh const &mesh, std::vector< MeshLib::Element * > &vec_matrix_elements, std::vector< int > &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Element * > > &vec_fracture_elements, std::vector< std::vector< MeshLib::Element * > > &vec_fracture_matrix_elements, std::vector< std::vector< MeshLib::Node * > > &vec_fracture_nodes, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_junction_nodeID_matIDs)
std::vector< double > uGlobalEnrichments(std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
BranchProperty createBranchProperty(MeshLib::Node const &branchNode, FractureProperty const &master_frac, FractureProperty const &slave_frac)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)