OGS
SmallDeformationProcess.cpp
Go to the documentation of this file.
1
12
14#include "MeshLib/Mesh.h"
15#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 std::vector<std::pair<std::size_t, std::vector<int>>>
48 vec_branch_nodeID_matIDs;
49 std::vector<std::pair<std::size_t, std::vector<int>>>
50 vec_junction_nodeID_matIDs;
54 _vec_fracture_nodes, vec_branch_nodeID_matIDs,
55 vec_junction_nodeID_matIDs);
56
57 if (_vec_fracture_mat_IDs.size() !=
58 _process_data.fracture_properties.size())
59 {
61 "The number of the given fracture properties ({:d}) are not "
62 "consistent with the number of fracture groups in a mesh ({:d}).",
63 _process_data.fracture_properties.size(),
65 }
66
67 // create a map from a material ID to a fracture ID
68 auto max_frac_mat_id = std::max_element(_vec_fracture_mat_IDs.begin(),
70 _process_data._map_materialID_to_fractureID.resize(*max_frac_mat_id + 1);
71 for (unsigned i = 0; i < _vec_fracture_mat_IDs.size(); i++)
72 {
73 _process_data._map_materialID_to_fractureID[_vec_fracture_mat_IDs[i]] =
74 i;
75 }
76
77 // create a table of connected fracture IDs for each element
78 _process_data._vec_ele_connected_fractureIDs.resize(
79 mesh.getNumberOfElements());
80 for (unsigned i = 0; i < _vec_fracture_matrix_elements.size(); i++)
81 {
82 for (auto e : _vec_fracture_matrix_elements[i])
83 {
84 _process_data._vec_ele_connected_fractureIDs[e->getID()].push_back(
85 i);
86 }
87 }
88
89 // set fracture property
90 for (auto& fracture_prop : _process_data.fracture_properties)
91 {
92 // based on the 1st element assuming a fracture forms a straight line
94 DisplacementDim,
95 *_vec_fracture_elements[fracture_prop.fracture_id][0],
96 fracture_prop);
97 }
98
99 // set branches
100 for (auto const& [vec_branch_nodeID, matID] : vec_branch_nodeID_matIDs)
101 {
102 auto master_matId = matID[0];
103 auto slave_matId = matID[1];
104 auto& master_frac =
105 _process_data.fracture_properties
106 [_process_data._map_materialID_to_fractureID[master_matId]];
107 auto& slave_frac =
108 _process_data.fracture_properties
109 [_process_data._map_materialID_to_fractureID[slave_matId]];
110
111 master_frac.branches_master.push_back(createBranchProperty(
112 *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
113
114 slave_frac.branches_slave.push_back(createBranchProperty(
115 *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
116 }
117
118 // set junctions
119 transform(cbegin(vec_junction_nodeID_matIDs),
120 cend(vec_junction_nodeID_matIDs),
121 back_inserter(_vec_junction_nodes),
122 [&](auto& vec_junction_nodeID_matID)
123 {
124 return const_cast<MeshLib::Node*>(
125 _mesh.getNode(vec_junction_nodeID_matID.first));
126 });
127
128 for (std::size_t i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
129 {
130 auto const& material_ids = vec_junction_nodeID_matIDs[i].second;
131 assert(material_ids.size() == 2);
132 std::array<int, 2> fracture_ids{
133 {_process_data._map_materialID_to_fractureID[material_ids[0]],
134 _process_data._map_materialID_to_fractureID[material_ids[1]]}};
135
136 _process_data.junction_properties.emplace_back(
137 i, *mesh.getNode(vec_junction_nodeID_matIDs[i].first),
138 fracture_ids);
139 }
140
141 // create a table of connected junction IDs for each element
142 _process_data._vec_ele_connected_junctionIDs.resize(
143 mesh.getNumberOfElements());
144 for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
145 {
146 auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
147 for (auto id :
149 {
150 _process_data._vec_ele_connected_junctionIDs[id].push_back(i);
151 }
152 }
153
154 // create a table of junction node and connected elements
156 vec_junction_nodeID_matIDs.size());
157 for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
158 {
159 auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
160 for (auto e : mesh.getElementsConnectedToNode(*node))
161 {
163 const_cast<MeshLib::Element*>(e));
164 }
165 }
166
167 //
168 // If Neumann BCs for the displacement_jump variable are required they need
169 // special treatment because of the levelset function. The implementation
170 // exists in the version 6.1.0 (e54815cc07ee89c81f953a4955b1c788595dd725)
171 // and was removed due to lack of applications.
172 //
173
174 MeshLib::PropertyVector<int> const* material_ids(
175 mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
176 _process_data._mesh_prop_materialIDs = material_ids;
177}
178
179template <int DisplacementDim>
181{
182 //------------------------------------------------------------
183 // prepare mesh subsets to define DoFs
184 //------------------------------------------------------------
185 // for extrapolation
186 _mesh_subset_all_nodes =
187 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
188 // regular u
189 _mesh_subset_matrix_nodes =
190 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
191 // u jump
192 for (unsigned i = 0; i < _vec_fracture_nodes.size(); i++)
193 {
194 _mesh_subset_fracture_nodes.push_back(
195 std::make_unique<MeshLib::MeshSubset const>(
196 _mesh, _vec_fracture_nodes[i]));
197 }
198 // enrichment for junctions
199 _mesh_subset_junction_nodes =
200 std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_junction_nodes);
201
202 // Collect the mesh subsets in a vector.
203 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
204 std::generate_n(std::back_inserter(all_mesh_subsets), DisplacementDim,
205 [&]() { return *_mesh_subset_matrix_nodes; });
206 for (auto const& ms : _mesh_subset_fracture_nodes)
207 {
208 std::generate_n(std::back_inserter(all_mesh_subsets),
209 DisplacementDim,
210 [&]() { return *ms; });
211 }
212 std::generate_n(std::back_inserter(all_mesh_subsets),
213 DisplacementDim,
214 [&]() { return *_mesh_subset_junction_nodes; });
215
216 std::vector<int> const vec_n_components(
217 1 + _vec_fracture_mat_IDs.size() + _vec_junction_nodes.size(),
218 DisplacementDim);
219
220 std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
221 vec_var_elements.push_back(&_vec_matrix_elements);
222 for (unsigned i = 0; i < _vec_fracture_matrix_elements.size(); i++)
223 {
224 vec_var_elements.push_back(&_vec_fracture_matrix_elements[i]);
225 }
226 for (unsigned i = 0; i < _vec_junction_fracture_matrix_elements.size(); i++)
227 {
228 vec_var_elements.push_back(&_vec_junction_fracture_matrix_elements[i]);
229 }
230
231 _local_to_global_index_map =
232 std::make_unique<NumLib::LocalToGlobalIndexMap>(
233 std::move(all_mesh_subsets),
234 vec_n_components,
235 vec_var_elements,
237}
238
239template <int DisplacementDim>
241 NumLib::LocalToGlobalIndexMap const& dof_table,
242 MeshLib::Mesh const& mesh,
243 unsigned const integration_order)
244{
249 mesh.getElements(), dof_table, _local_assemblers,
250 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
251 _process_data);
252
253 // TODO move the two data members somewhere else.
254 // for extrapolation of secondary variables
255 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
256 *_mesh_subset_all_nodes};
257 _local_to_global_index_map_single_component =
258 std::make_unique<NumLib::LocalToGlobalIndexMap>(
259 std::move(all_mesh_subsets_single_component),
260 // by location order is needed for output
262
263 auto add_secondary_variable = [&](std::string const& name,
264 int const num_components,
265 auto get_ip_values_function)
266 {
267 _secondary_variables.addSecondaryVariable(
268 name,
269 makeExtrapolator(num_components, getExtrapolator(),
270 _local_assemblers,
271 std::move(get_ip_values_function)));
272 };
273
274 add_secondary_variable("sigma",
276 DisplacementDim>::RowsAtCompileTime,
277 &LocalAssemblerInterface::getIntPtSigma);
278
279 add_secondary_variable("epsilon",
281 DisplacementDim>::RowsAtCompileTime,
282 &LocalAssemblerInterface::getIntPtEpsilon);
283
284 add_secondary_variable("fracture_stress", DisplacementDim,
285 &LocalAssemblerInterface::getIntPtFractureStress);
286
287 add_secondary_variable("fracture_aperture", 1,
288 &LocalAssemblerInterface::getIntPtFractureAperture);
289
290 _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
291 const_cast<MeshLib::Mesh&>(mesh), "sigma_avg",
294 DisplacementDim>::RowsAtCompileTime);
295
296 for (MeshLib::Element const* e : _mesh.getElements())
297 {
298 if (e->getDimension() < DisplacementDim)
299 {
300 continue;
301 }
302
303 Eigen::Vector3d const pt(getCenterOfGravity(*e).asEigenVector3d());
304 std::vector<FractureProperty*> e_fracture_props;
305 std::unordered_map<int, int> e_fracID_to_local;
306 unsigned tmpi = 0;
307 for (auto fid :
308 _process_data._vec_ele_connected_fractureIDs[e->getID()])
309 {
310 e_fracture_props.push_back(&_process_data.fracture_properties[fid]);
311 e_fracID_to_local.insert({fid, tmpi++});
312 }
313 std::vector<JunctionProperty*> e_junction_props;
314 std::unordered_map<int, int> e_juncID_to_local;
315 tmpi = 0;
316 for (auto fid :
317 _process_data._vec_ele_connected_junctionIDs[e->getID()])
318 {
319 e_junction_props.push_back(&_process_data.junction_properties[fid]);
320 e_juncID_to_local.insert({fid, tmpi++});
321 }
322 std::vector<double> const levelsets(uGlobalEnrichments(
323 e_fracture_props, e_junction_props, e_fracID_to_local, pt));
324
325 for (unsigned i = 0; i < e_fracture_props.size(); i++)
326 {
327 auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
328 const_cast<MeshLib::Mesh&>(mesh),
329 "levelset" +
330 std::to_string(e_fracture_props[i]->fracture_id + 1),
332 mesh_prop_levelset->resize(mesh.getNumberOfElements());
333 (*mesh_prop_levelset)[e->getID()] = levelsets[i];
334 }
335 for (unsigned i = 0; i < e_junction_props.size(); i++)
336 {
337 auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
338 const_cast<MeshLib::Mesh&>(mesh),
339 "levelset" +
340 std::to_string(e_junction_props[i]->junction_id + 1 +
341 _process_data.fracture_properties.size()),
343 mesh_prop_levelset->resize(mesh.getNumberOfElements());
344 (*mesh_prop_levelset)[e->getID()] =
345 levelsets[i + e_fracture_props.size()];
346 }
347 }
348
349 _process_data.element_local_jumps =
350 MeshLib::getOrCreateMeshProperty<double>(
351 const_cast<MeshLib::Mesh&>(mesh), "local_jump_w_avg",
352 MeshLib::MeshItemType::Cell, DisplacementDim);
353
354 _process_data.element_fracture_stresses =
355 MeshLib::getOrCreateMeshProperty<double>(
356 const_cast<MeshLib::Mesh&>(mesh), "fracture_stress_avg",
357 MeshLib::MeshItemType::Cell, DisplacementDim);
358
359 auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
360 const_cast<MeshLib::Mesh&>(mesh), "fracture_aperture_avg",
362
363 mesh_prop_b->resize(mesh.getNumberOfElements());
364 auto const& mesh_prop_matid = *_process_data._mesh_prop_materialIDs;
365 for (auto const& fracture_prop : _process_data.fracture_properties)
366 {
367 for (MeshLib::Element const* e : _mesh.getElements())
368 {
369 if (e->getDimension() == DisplacementDim)
370 {
371 continue;
372 }
373 if (mesh_prop_matid[e->getID()] != fracture_prop.mat_id)
374 {
375 continue;
376 }
377 // Mean value for the element. This allows usage of node based
378 // properties for aperture.
379 (*mesh_prop_b)[e->getID()] =
380 fracture_prop.aperture0
381 .getNodalValuesOnElement(*e, /*time independent*/ 0)
382 .mean();
383 }
384 }
385 _process_data._mesh_prop_b = mesh_prop_b;
386}
387
388template <int DisplacementDim>
390 double const t, double const dt, std::vector<GlobalVector*> const& x,
391 GlobalVector const& x_prev, int const process_id)
392{
393 DBUG("Compute the secondary variables for SmallDeformationProcess.");
394 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
395 dof_tables.reserve(x.size());
396 std::generate_n(std::back_inserter(dof_tables), x.size(),
397 [&]() { return _local_to_global_index_map.get(); });
398
399 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
402 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
403 x_prev, process_id);
404}
405
406template <int DisplacementDim>
408{
409 return false;
410}
411
412template <int DisplacementDim>
414 const double t, double const dt, std::vector<GlobalVector*> const& x,
415 std::vector<GlobalVector*> const& x_prev, int const process_id,
417{
418 DBUG("Assemble SmallDeformationProcess.");
419
420 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
421
422 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
423 dof_table = {std::ref(*_local_to_global_index_map)};
424 // Call global assembler for each local assembly item.
426 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
427 pv.getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, M, K,
428 b);
429}
430template <int DisplacementDim>
433 const double t, double const dt, std::vector<GlobalVector*> const& x,
434 std::vector<GlobalVector*> const& x_prev, int const process_id,
436{
437 DBUG("AssembleWithJacobian SmallDeformationProcess.");
438
439 // Call global assembler for each local assembly item.
440 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
441
442 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
443 dof_table = {std::ref(*_local_to_global_index_map)};
446 _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x,
447 x_prev, process_id, M, K, b, Jac);
448}
449template <int DisplacementDim>
451 std::vector<GlobalVector*> const& x, double const t, double const dt,
452 const int process_id)
453{
454 DBUG("PreTimestep SmallDeformationProcess.");
455
456 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
459 _local_assemblers, pv.getActiveElementIDs(),
460 *_local_to_global_index_map, *x[process_id], t, dt);
461}
462// ------------------------------------------------------------------------------------
463// template instantiation
464// ------------------------------------------------------------------------------------
465template class SmallDeformationProcess<2>;
466template class SmallDeformationProcess<3>;
467
468} // namespace SmallDeformation
469} // namespace LIE
470} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
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:25
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:63
bool isAxiallySymmetric() const
Definition Mesh.h:137
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition Mesh.h:91
Properties & getProperties()
Definition Mesh.h:134
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:256
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:97
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 assembleWithJacobianConcreteProcess(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, GlobalMatrix &Jac) override
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
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 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)
std::vector< std::size_t > const & getActiveElementIDs() const
MeshLib::Mesh & _mesh
Definition Process.h:354
Handles configuration of several secondary variables from the project file.
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, 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, GlobalMatrix &Jac)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > 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)
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:225
@ 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)