Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "MeshGeneratorPD.h"
11 : #include "PeridynamicsMesh.h"
12 : #include "CastUniquePointer.h"
13 :
14 : #include "libmesh/edge_edge2.h"
15 : #include "libmesh/face_tri3.h"
16 : #include "libmesh/cell_tet4.h"
17 :
18 : registerMooseObject("PeridynamicsApp", MeshGeneratorPD);
19 :
20 : InputParameters
21 632 : MeshGeneratorPD::validParams()
22 : {
23 632 : InputParameters params = MeshGenerator::validParams();
24 632 : params.addClassDescription("Mesh generator class to convert FE mesh to PD mesh");
25 :
26 1264 : params.addRequiredParam<MeshGeneratorName>("input",
27 : "The mesh based on which PD mesh will be created");
28 1264 : params.addParam<std::vector<SubdomainID>>("blocks_to_pd",
29 : "IDs of the FE mesh blocks to be converted to PD mesh");
30 1264 : params.addParam<std::vector<SubdomainID>>(
31 : "blocks_as_fe",
32 : "IDs of the FE mesh blocks to not be converted to PD mesh. This should only be used when the "
33 : "number of to-be-converted FE blocks is considerably large.");
34 1264 : params.addRequiredParam<bool>(
35 : "retain_fe_mesh", "Whether to retain the FE mesh or not after conversion into PD mesh");
36 1264 : params.addParam<bool>(
37 : "construct_pd_sidesets",
38 1264 : false,
39 : "Whether to construct PD sidesets based on the sidesets in original FE mesh");
40 1264 : params.addParam<std::vector<boundary_id_type>>(
41 : "sidesets_to_pd", "IDs of the FE sidesets to be reconstructed based on converted PD mesh");
42 1264 : params.addParam<std::vector<std::vector<SubdomainID>>>(
43 : "bonding_block_pairs",
44 : "List of FE block pairs between which inter-block bonds will be created after being "
45 : "converted into PD mesh");
46 1264 : params.addParam<std::vector<std::vector<SubdomainID>>>(
47 : "non_bonding_block_pairs",
48 : "List of FE block pairs between which inter-block bonds will NOT be created after being "
49 : "converted into PD mesh");
50 1264 : params.addParam<bool>(
51 : "merge_pd_blocks",
52 1264 : false,
53 : "Whether to merge all converted PD mesh blocks into a single block. This is "
54 : "used when all PD blocks have the same properties");
55 1264 : params.addParam<bool>(
56 : "merge_pd_interfacial_blocks",
57 1264 : false,
58 : "Whether to merge all PD interfacial mesh blocks into a single block. This is used "
59 : "when all PD interfacial blocks have the same properties");
60 :
61 632 : return params;
62 0 : }
63 :
64 316 : MeshGeneratorPD::MeshGeneratorPD(const InputParameters & parameters)
65 : : MeshGenerator(parameters),
66 316 : _input(getMesh("input")),
67 632 : _has_blks_to_pd(isParamValid("blocks_to_pd")),
68 632 : _has_blks_as_fe(isParamValid("blocks_as_fe")),
69 632 : _retain_fe_mesh(getParam<bool>("retain_fe_mesh")),
70 632 : _merge_pd_blks(getParam<bool>("merge_pd_blocks")),
71 632 : _construct_pd_sideset(getParam<bool>("construct_pd_sidesets")),
72 632 : _has_sidesets_to_pd(isParamValid("sidesets_to_pd")),
73 632 : _has_bonding_blk_pairs(isParamValid("bonding_block_pairs")),
74 632 : _has_non_bonding_blk_pairs(isParamValid("non_bonding_block_pairs")),
75 948 : _merge_pd_interfacial_blks(getParam<bool>("merge_pd_interfacial_blocks"))
76 : {
77 316 : if (_has_blks_to_pd && _has_blks_as_fe)
78 0 : mooseError("Please specifiy either 'blocks_to_pd' or 'blocks_as_fe'!");
79 :
80 316 : if (_has_blks_as_fe)
81 : {
82 0 : std::vector<SubdomainID> ids = getParam<std::vector<SubdomainID>>("blocks_as_fe");
83 0 : for (unsigned int i = 0; i < ids.size(); ++i)
84 0 : _blks_as_fe.insert(ids[i]);
85 0 : }
86 :
87 316 : if (_has_blks_to_pd)
88 : {
89 44 : std::vector<SubdomainID> ids = getParam<std::vector<SubdomainID>>("blocks_to_pd");
90 64 : for (unsigned int i = 0; i < ids.size(); ++i)
91 42 : _blks_to_pd.insert(ids[i]);
92 22 : }
93 :
94 316 : if (!_construct_pd_sideset && _has_sidesets_to_pd)
95 0 : mooseError("'sidesets_to_pd' is provided without setting "
96 : "'construct_pd_sidesets' to 'true'!");
97 :
98 316 : if (_has_sidesets_to_pd)
99 : {
100 4 : std::vector<boundary_id_type> ids = getParam<std::vector<boundary_id_type>>("sidesets_to_pd");
101 4 : for (unsigned int i = 0; i < ids.size(); ++i)
102 2 : _fe_sidesets_for_pd_construction.insert(ids[i]);
103 2 : }
104 :
105 316 : if (_has_bonding_blk_pairs && _has_non_bonding_blk_pairs)
106 0 : mooseError("Please specifiy either 'bonding_block_pairs' or "
107 : "'non_bonding_block_pairs'!");
108 :
109 : _pd_bonding_blk_pairs.clear();
110 316 : if (_has_bonding_blk_pairs)
111 : {
112 : std::vector<std::vector<SubdomainID>> id_pairs =
113 14 : getParam<std::vector<std::vector<SubdomainID>>>("bonding_block_pairs");
114 :
115 24 : for (unsigned int i = 0; i < id_pairs.size(); ++i)
116 17 : _pd_bonding_blk_pairs.insert(std::make_pair(id_pairs[i][0] + _pd_blk_offset_number,
117 17 : id_pairs[i][1] + _pd_blk_offset_number));
118 7 : } // considered the renumbering of IDs of converted blocks
119 :
120 : _pd_non_bonding_blk_pairs.clear();
121 316 : if (_has_non_bonding_blk_pairs)
122 : {
123 : std::vector<std::vector<SubdomainID>> id_pairs =
124 0 : getParam<std::vector<std::vector<SubdomainID>>>("non_bonding_block_pairs");
125 :
126 0 : for (unsigned int i = 0; i < id_pairs.size(); ++i)
127 0 : _pd_non_bonding_blk_pairs.insert(std::make_pair(id_pairs[i][0] + _pd_blk_offset_number,
128 0 : id_pairs[i][1] + _pd_blk_offset_number));
129 0 : } // considered the renumbering of IDs of converted blocks
130 316 : }
131 :
132 : std::unique_ptr<MeshBase>
133 316 : MeshGeneratorPD::generate()
134 : {
135 : // get the MeshBase object this generator will be working on
136 316 : std::unique_ptr<MeshBase> old_mesh = std::move(_input);
137 :
138 : // If it's not prepared, we need it to be, to access caches like
139 : // mesh_dimension() and spatial_dimension()
140 316 : if (!old_mesh->is_prepared())
141 250 : old_mesh->prepare_for_use();
142 :
143 : // STEP 1: obtain FE block(s) and elements to be converted to PD mesh
144 :
145 : // get the IDs of all available blocks in the input FE mesh
146 : std::set<SubdomainID> all_fe_blks;
147 100916 : for (const auto & old_elem : old_mesh->element_ptr_range())
148 50458 : all_fe_blks.insert(old_elem->subdomain_id());
149 :
150 : // double check the existence of the block ids provided by user
151 316 : if (_has_blks_to_pd)
152 62 : for (const auto & blkit : _blks_to_pd)
153 : if (!all_fe_blks.count(blkit))
154 2 : mooseError("Block ID ", blkit, " in the 'blocks_to_pd' does not exist in the FE mesh!");
155 :
156 314 : if (_has_blks_as_fe)
157 0 : for (const auto & blkit : _blks_as_fe)
158 : if (!all_fe_blks.count(blkit))
159 0 : mooseError("Block ID ", blkit, " in the 'blocks_as_fe' does not exist in the FE mesh!");
160 :
161 : // the maximum FE block ID, which will be used in determine the block ID for interfacial bond in
162 : // the case of single interface block
163 314 : const unsigned int max_fe_blk_id = *all_fe_blks.rbegin();
164 :
165 : // categorize mesh blocks into converted and non-converted blocks
166 314 : if (_has_blks_to_pd)
167 20 : std::set_difference(all_fe_blks.begin(),
168 : all_fe_blks.end(),
169 : _blks_to_pd.begin(),
170 : _blks_to_pd.end(),
171 20 : std::inserter(_blks_as_fe, _blks_as_fe.begin()));
172 294 : else if (_has_blks_as_fe)
173 0 : std::set_difference(all_fe_blks.begin(),
174 : all_fe_blks.end(),
175 : _blks_as_fe.begin(),
176 : _blks_as_fe.end(),
177 0 : std::inserter(_blks_to_pd, _blks_to_pd.begin()));
178 : else // if no block ids provided by user, by default, convert all FE mesh to PD mesh
179 : _blks_to_pd = all_fe_blks;
180 :
181 314 : if (_has_bonding_blk_pairs)
182 : {
183 22 : for (const auto & blk_pair : _pd_bonding_blk_pairs)
184 : {
185 17 : if (!all_fe_blks.count(blk_pair.first - _pd_blk_offset_number))
186 2 : mooseError("Block ID ",
187 2 : blk_pair.first - _pd_blk_offset_number,
188 : " in the 'bonding_block_pairs' does not exist in the FE mesh!");
189 15 : if (!all_fe_blks.count(blk_pair.second - _pd_blk_offset_number))
190 0 : mooseError("Block ID ",
191 0 : blk_pair.second - _pd_blk_offset_number,
192 : " in the 'bonding_block_pairs' does not exist in the FE mesh!");
193 15 : if (!_blks_to_pd.count(blk_pair.first - _pd_blk_offset_number))
194 0 : mooseError("Block ID ",
195 0 : blk_pair.first - _pd_blk_offset_number,
196 : " in the 'bonding_block_pairs' is a FE mesh block!");
197 15 : if (!_blks_to_pd.count(blk_pair.second - _pd_blk_offset_number))
198 0 : mooseError("Block ID ",
199 0 : blk_pair.second - _pd_blk_offset_number,
200 : " in the 'bonding_block_pairs' is a FE mesh block!");
201 : }
202 : }
203 :
204 312 : if (_has_non_bonding_blk_pairs)
205 0 : for (const auto & blk_pair : _pd_non_bonding_blk_pairs)
206 : {
207 0 : if (!all_fe_blks.count(blk_pair.first - _pd_blk_offset_number))
208 0 : mooseError("Block ID ",
209 0 : blk_pair.first - _pd_blk_offset_number,
210 : " in the 'non_bonding_block_pairs' does not exist in the FE mesh!");
211 0 : if (!all_fe_blks.count(blk_pair.second - _pd_blk_offset_number))
212 0 : mooseError("Block ID ",
213 0 : blk_pair.second - _pd_blk_offset_number,
214 : " in the 'non_bonding_block_pairs' does not exist in the FE mesh!");
215 0 : if (!_blks_as_fe.count(blk_pair.first - _pd_blk_offset_number))
216 0 : mooseError("Block ID ",
217 0 : blk_pair.first - _pd_blk_offset_number,
218 : " in the 'non_bonding_block_pairs' is a FE mesh block!");
219 0 : if (!_blks_as_fe.count(blk_pair.second - _pd_blk_offset_number))
220 0 : mooseError("Block ID ",
221 0 : blk_pair.second - _pd_blk_offset_number,
222 : " in the 'non_bonding_block_pairs' is a FE mesh block!");
223 : }
224 :
225 : // the minimum converted FE block ID, which will be used to assign block ID for non-interfacial
226 : // bond in the case of combine converted blocks
227 312 : const unsigned int min_converted_fe_blk_id = *_blks_to_pd.begin();
228 :
229 : // IDs of to-be-converted FE elems
230 : std::set<dof_id_type> elems_to_pd;
231 : // retained FE mesh and non-converted FE mesh, if any
232 : std::set<dof_id_type> fe_nodes;
233 : std::set<dof_id_type> fe_elems;
234 100060 : for (const auto & old_elem : old_mesh->element_ptr_range())
235 49718 : if (_blks_to_pd.count(old_elem->subdomain_id())) // record to-be-converted FE elem IDs
236 : {
237 49188 : elems_to_pd.insert(old_elem->id());
238 49188 : if (_retain_fe_mesh) // save converted elems and their nodes if need to be retained
239 : {
240 730 : fe_elems.insert(old_elem->id());
241 2920 : for (unsigned int i = 0; i < old_elem->n_nodes(); ++i)
242 2190 : fe_nodes.insert(old_elem->node_id(i));
243 : }
244 : }
245 : else // save non-converted elements and their nodes
246 : {
247 530 : fe_elems.insert(old_elem->id());
248 2120 : for (unsigned int i = 0; i < old_elem->n_nodes(); ++i)
249 1590 : fe_nodes.insert(old_elem->node_id(i));
250 312 : }
251 :
252 : // number of FE elements and nodes in the old mesh to be saved in the new mesh
253 : dof_id_type n_fe_nodes = fe_nodes.size();
254 : dof_id_type n_fe_elems = fe_elems.size();
255 : dof_id_type n_phantom_elems = 0;
256 :
257 : // determine the number of phantom elements to be generated in the new mesh based on sideset in
258 : // old mesh
259 : BoundaryInfo & old_boundary_info = old_mesh->get_boundary_info();
260 : const std::set<boundary_id_type> & all_fe_sidesets = old_boundary_info.get_side_boundary_ids();
261 :
262 : // double check the existence of the sideset ids provided by user
263 312 : if (_has_sidesets_to_pd)
264 : {
265 2 : for (const auto & sideset : _fe_sidesets_for_pd_construction)
266 : if (!all_fe_sidesets.count(sideset))
267 2 : mooseError("Sideset ID ",
268 : sideset,
269 : " in the 'sidesets_to_pd' does not exist in the finite "
270 : "element mesh!");
271 : }
272 : else // save the IDs of all FE sidesets, this will be updated to the converted FE block later
273 : _fe_sidesets_for_pd_construction = all_fe_sidesets;
274 :
275 : // determine number of FE side elements, the number of actual phantom elements is less than or
276 : // equal to the number of FE side elements, this number is used to reserve number of elements
277 : // in the new mesh only
278 : std::map<boundary_id_type, std::set<dof_id_type>> fe_bid_eid;
279 310 : auto fe_sbc_tuples = old_boundary_info.build_side_list();
280 : // 0: element ID, 1: side ID, 2: boundary ID
281 9655 : for (const auto & sbct : fe_sbc_tuples)
282 : if (_fe_sidesets_for_pd_construction.count(
283 : std::get<2>(sbct))) // save elements of to-be-constructed sidesets only
284 : {
285 9345 : fe_bid_eid[std::get<2>(sbct)].insert(std::get<0>(sbct));
286 9345 : ++n_phantom_elems;
287 : }
288 :
289 : // STEP 2: generate PD data based on to-be-converted FE mesh and prepare for new mesh
290 :
291 310 : PeridynamicsMesh & pd_mesh = dynamic_cast<PeridynamicsMesh &>(*_mesh);
292 : // generate PD node data, including neighbors
293 620 : pd_mesh.createPeridynamicsMeshData(
294 : *old_mesh, elems_to_pd, _pd_bonding_blk_pairs, _pd_non_bonding_blk_pairs);
295 :
296 : // number of PD elements and nodes to be created
297 310 : dof_id_type n_pd_nodes = pd_mesh.nPDNodes();
298 310 : dof_id_type n_pd_bonds = pd_mesh.nPDBonds();
299 :
300 : // initialize an empty new mesh object
301 310 : std::unique_ptr<MeshBase> new_mesh = buildMeshBaseObject();
302 310 : new_mesh->clear();
303 : // set new mesh dimension
304 620 : new_mesh->set_mesh_dimension(old_mesh->mesh_dimension());
305 310 : new_mesh->set_spatial_dimension(old_mesh->spatial_dimension());
306 : // reserve elements and nodes for the new mesh
307 310 : new_mesh->reserve_nodes(n_pd_nodes + n_fe_nodes);
308 310 : new_mesh->reserve_elem(n_pd_bonds + n_fe_elems + n_phantom_elems);
309 :
310 : BoundaryInfo & new_boundary_info = new_mesh->get_boundary_info();
311 :
312 : // STEP 3: add points of PD and FE (retained and/or non-converted) nodes, if any, to new mesh
313 :
314 : // save PD nodes to new mesh first
315 : unsigned int new_node_id = 0;
316 : // map of IDs of converted FE elements and PD nodes
317 : std::map<dof_id_type, dof_id_type> fe_elem_pd_node_map;
318 49110 : for (const auto & eid : elems_to_pd)
319 : {
320 48800 : new_mesh->add_point(old_mesh->elem_ptr(eid)->vertex_average(), new_node_id);
321 48800 : fe_elem_pd_node_map.insert(std::make_pair(eid, new_node_id));
322 :
323 48800 : ++new_node_id;
324 : }
325 : // then save both retained and non-converted FE nodes, if any, to the new mesh
326 : // map of IDs of the same point in old and new meshes
327 : std::map<dof_id_type, dof_id_type> fe_nodes_map;
328 1350 : for (const auto & nid : fe_nodes)
329 : {
330 1040 : new_mesh->add_point(*old_mesh->node_ptr(nid), new_node_id);
331 1040 : fe_nodes_map.insert(std::make_pair(old_mesh->node_ptr(nid)->id(), new_node_id));
332 :
333 1040 : ++new_node_id;
334 : }
335 :
336 : // STEP 4: generate PD, phantom, and FE elems using retained and/or non-converted meshes if any
337 :
338 : // first, generate PD elements for new mesh
339 : unsigned int new_elem_id = 0;
340 49110 : for (unsigned int i = 0; i < n_pd_nodes; ++i)
341 : {
342 48800 : std::vector<dof_id_type> pd_node_neighbors = pd_mesh.getNeighbors(i); // neighbors of a PD node
343 2094570 : for (unsigned int j = 0; j < pd_node_neighbors.size(); ++j)
344 2045770 : if (pd_node_neighbors[j] > i)
345 : {
346 1022885 : SubdomainID bid_i = pd_mesh.getNodeBlockID(i);
347 1022885 : SubdomainID bid_j = pd_mesh.getNodeBlockID(pd_node_neighbors[j]);
348 1022885 : Elem * new_elem = new Edge2;
349 1022885 : new_elem->set_id(new_elem_id);
350 1022885 : if (bid_i == bid_j) // assign block ID to PD non-inter-block bonds
351 1017790 : if (_merge_pd_blks)
352 33865 : new_elem->subdomain_id() = min_converted_fe_blk_id + _pd_blk_offset_number;
353 : else
354 983925 : new_elem->subdomain_id() = bid_i;
355 5095 : else if (_merge_pd_interfacial_blks) // assign block ID (max_fe_blk_id + 1 +
356 : // _pd_blk_offset_number) to all PD inter-block bonds
357 5095 : new_elem->subdomain_id() = max_fe_blk_id + 1 + _pd_blk_offset_number;
358 : else // assign a new block ID (node i blk ID + node j blk ID) to this PD inter-block bonds
359 0 : new_elem->subdomain_id() = bid_i + bid_j;
360 :
361 1022885 : new_elem = new_mesh->add_elem(new_elem);
362 1022885 : new_elem->set_node(0, new_mesh->node_ptr(i));
363 1022885 : new_elem->set_node(1, new_mesh->node_ptr(pd_node_neighbors[j]));
364 :
365 1022885 : ++new_elem_id;
366 : }
367 48800 : }
368 :
369 310 : if (_merge_pd_blks) // update PD node block ID if use single block ID for all PD blocks
370 5 : pd_mesh.setNodeBlockID(min_converted_fe_blk_id + _pd_blk_offset_number);
371 :
372 : // then generate phantom elements for sidesets in PD mesh, this is optional
373 : std::map<std::pair<dof_id_type, dof_id_type>, std::set<dof_id_type>> elem_edge_nodes;
374 310 : if (_construct_pd_sideset)
375 : {
376 75 : for (const auto & bidit : _fe_sidesets_for_pd_construction)
377 : {
378 1945 : for (const auto & eidit : fe_bid_eid[bidit])
379 : {
380 : bool should_add = false;
381 1890 : Elem * old_elem = old_mesh->elem_ptr(eidit);
382 : if (_blks_to_pd.count(
383 : old_elem->subdomain_id())) // limit the sidesets to the converted FE blocks only
384 : {
385 : std::vector<dof_id_type> pd_nodes;
386 : Point p0, p1, p2, p3;
387 :
388 1890 : if (old_elem->dim() == 2) // 2D: construct 3-node triangular phanton elems
389 : {
390 330 : pd_nodes.resize(3);
391 : // the current elem on the boundary side is the first node of a phantom element
392 330 : pd_nodes[0] = fe_elem_pd_node_map.at(eidit);
393 330 : p0 = *new_mesh->node_ptr(pd_nodes[0]);
394 :
395 : // search for two more nodes to construct a phantom 3-node triangular element
396 430 : if (old_elem->n_nodes() == 3 ||
397 100 : old_elem->n_nodes() == 6) // original triangular FE elements: 3-node or 6-node
398 : {
399 : // one special case of triangular mesh: two elems share a boundary node
400 : bool has_neighbor_on_boundary = false;
401 920 : for (unsigned int i = 0; i < old_elem->n_neighbors(); ++i)
402 : {
403 : Elem * nb = old_elem->neighbor_ptr(i);
404 690 : if (nb != nullptr)
405 : {
406 450 : if (fe_bid_eid[bidit].count(nb->id()))
407 : {
408 : has_neighbor_on_boundary = true;
409 30 : pd_nodes[1] = fe_elem_pd_node_map.at(nb->id());
410 30 : p1 = *new_mesh->node_ptr(pd_nodes[1]);
411 : }
412 : else
413 : {
414 420 : pd_nodes[2] = fe_elem_pd_node_map.at(nb->id());
415 420 : p2 = *new_mesh->node_ptr(pd_nodes[2]);
416 : }
417 : }
418 : }
419 :
420 230 : if (has_neighbor_on_boundary &&
421 30 : (p1(1) - p0(1)) * (p2(0) - p1(0)) - (p1(0) - p0(0)) * (p2(1) - p1(1)) < 0)
422 : {
423 : should_add = true;
424 : }
425 : else // common case of triangular mesh: three elems share a boundary node
426 : {
427 860 : for (unsigned int i = 0; i < old_elem->n_neighbors(); ++i)
428 : {
429 : Elem * nb = old_elem->neighbor_ptr(i);
430 :
431 645 : if (nb != nullptr)
432 : {
433 420 : p2 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(nb->id()));
434 :
435 2100 : for (unsigned int j = 0; j < nb->n_neighbors(); ++j)
436 : {
437 : Elem * nb_nb = nb->neighbor_ptr(j);
438 :
439 2475 : if (nb_nb != nullptr && fe_bid_eid[bidit].count(nb_nb->id()) &&
440 750 : nb_nb->id() != eidit)
441 : {
442 330 : p1 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(nb_nb->id()));
443 :
444 : // add new phantom elem only if (p0, p1, p2) is counterclockwisely ordered
445 330 : if ((p1(1) - p0(1)) * (p2(0) - p1(0)) - (p1(0) - p0(0)) * (p2(1) - p1(1)) <
446 : 0)
447 : {
448 : should_add = true;
449 170 : pd_nodes[1] = fe_elem_pd_node_map.at(nb_nb->id());
450 170 : pd_nodes[2] = fe_elem_pd_node_map.at(nb->id());
451 : }
452 : }
453 : }
454 :
455 420 : if (!should_add) // if using the neighbor's neighbor is not a success
456 : {
457 : // one special case of triangular mesh: four elems share a boundary node
458 : // other cases of more than four elems share a boundary node is not considered
459 300 : for (unsigned int i = 0; i < old_elem->n_neighbors(); ++i)
460 : {
461 : Elem * nb = old_elem->neighbor_ptr(i);
462 :
463 225 : if (nb != nullptr)
464 : {
465 145 : p2 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(nb->id()));
466 :
467 725 : for (unsigned int j = 0; j < nb->n_neighbors(); ++j)
468 : {
469 : Elem * nb_nb = nb->neighbor_ptr(j);
470 435 : if (nb_nb != nullptr && nb_nb->id() != eidit)
471 : {
472 1040 : for (unsigned int k = 0; k < nb_nb->n_neighbors(); ++k)
473 : {
474 : Elem * nb_nb_nb = nb_nb->neighbor_ptr(k);
475 :
476 : if (nb_nb_nb != nullptr &&
477 1470 : fe_bid_eid[bidit].count(nb_nb_nb->id()) &&
478 15 : nb_nb_nb->id() != eidit)
479 : {
480 15 : p1 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(nb_nb_nb->id()));
481 :
482 : // add new phantom elem only if (p0, p1, p2) is counterclockwisely
483 : // ordered
484 15 : if ((p1(1) - p0(1)) * (p2(0) - p1(0)) -
485 15 : (p1(0) - p0(0)) * (p2(1) - p1(1)) <
486 : 0)
487 : {
488 : should_add = true;
489 15 : pd_nodes[1] = fe_elem_pd_node_map.at(nb_nb_nb->id());
490 15 : pd_nodes[2] = fe_elem_pd_node_map.at(nb->id());
491 : }
492 : }
493 : }
494 : }
495 : }
496 : }
497 : }
498 : }
499 : }
500 : }
501 : }
502 : }
503 100 : else if (old_elem->n_nodes() == 4 || old_elem->n_nodes() == 8 ||
504 0 : old_elem->n_nodes() ==
505 : 9) // original quadrilateral FE elements: 4-node, 8-node and 9-node
506 : {
507 : std::vector<dof_id_type> pd_boundary_node_ids;
508 600 : for (unsigned int i = 0; i < old_elem->n_neighbors(); ++i)
509 : {
510 : Elem * nb = old_elem->neighbor_ptr(i);
511 500 : if (nb != nullptr && fe_bid_eid[bidit].count(nb->id()))
512 180 : pd_boundary_node_ids.push_back(fe_elem_pd_node_map.at(nb->id()));
513 : }
514 : // the neighbor opposite to the boundary side is the third node for the phantom elem
515 100 : pd_nodes[2] = fe_elem_pd_node_map.at(
516 : old_elem
517 100 : ->neighbor_ptr(old_elem->opposite_side(
518 100 : old_boundary_info.side_with_boundary_id(old_elem, bidit)))
519 100 : ->id());
520 100 : p2 = *new_mesh->node_ptr(pd_nodes[2]);
521 :
522 : // if two boundary neighbors, one of them is the second node for the phantom elem
523 : // if one boundary neighbors, it is the second node for the phantom elem
524 280 : for (unsigned int i = 0; i < pd_boundary_node_ids.size(); ++i)
525 : {
526 180 : p1 = *new_mesh->node_ptr(pd_boundary_node_ids[i]);
527 :
528 : // new phantom elem only if (p0, p1, p2) is counterclockwisely orderd
529 180 : if ((p1(1) - p0(1)) * (p2(0) - p1(0)) - (p1(0) - p0(0)) * (p2(1) - p1(1)) < 0)
530 : {
531 : should_add = true;
532 90 : pd_nodes[1] = pd_boundary_node_ids[i];
533 : }
534 : }
535 100 : }
536 : else
537 0 : mooseError("Element type ",
538 0 : old_elem->type(),
539 : " is not supported for PD sideset construction!");
540 :
541 315 : if (should_add)
542 : {
543 290 : Elem * new_elem = new Tri3;
544 290 : new_elem->set_id(new_elem_id);
545 290 : if (_merge_pd_blks)
546 220 : new_elem->subdomain_id() = min_converted_fe_blk_id + _phantom_blk_offset_number;
547 : else
548 70 : new_elem->subdomain_id() = old_elem->subdomain_id() + _phantom_blk_offset_number;
549 290 : new_elem = new_mesh->add_elem(new_elem);
550 290 : new_elem->set_node(0, new_mesh->node_ptr(pd_nodes[0]));
551 290 : new_elem->set_node(1, new_mesh->node_ptr(pd_nodes[1]));
552 290 : new_elem->set_node(2, new_mesh->node_ptr(pd_nodes[2]));
553 :
554 290 : ++new_elem_id;
555 :
556 290 : new_boundary_info.add_side(new_elem, 0, _pd_blk_offset_number + bidit);
557 290 : if (old_boundary_info.get_sideset_name(bidit) != "")
558 0 : new_boundary_info.sideset_name(_pd_blk_offset_number + bidit) =
559 0 : "pd_side_" + old_boundary_info.get_sideset_name(bidit);
560 : }
561 : }
562 1560 : else if (old_elem->dim() == 3) // 3D
563 : {
564 1560 : pd_nodes.resize(4); // construct 4-node tetrahedral phanton elems
565 1560 : pd_nodes[0] = fe_elem_pd_node_map.at(eidit);
566 1560 : p0 = *new_mesh->node_ptr(pd_nodes[0]);
567 : // search for three more nodes to construct a phantom 4-node tetrahedral element
568 3120 : if (old_elem->n_nodes() == 4 ||
569 1560 : old_elem->n_nodes() == 10) // original tetrahedral FE elements: 4-node or 10-node
570 : {
571 : // construct phantom element based on original tet mesh
572 0 : mooseError("Peridynamics sideset generation does not accept tetrahedral elements!");
573 : }
574 1560 : else if (old_elem->n_nodes() == 8 || old_elem->n_nodes() == 20 ||
575 0 : old_elem->n_nodes() ==
576 : 27) // original hexahedral FE elements: 8-node, 20-node or 27-node
577 : {
578 : std::vector<dof_id_type> boundary_elem_ids;
579 12480 : for (unsigned int i = 0; i < old_elem->n_neighbors(); ++i)
580 : {
581 : Elem * nb = old_elem->neighbor_ptr(i);
582 11280 : if (nb != nullptr && fe_bid_eid[bidit].count(nb->id()))
583 5520 : boundary_elem_ids.push_back(nb->id());
584 : }
585 :
586 : // find the fourth pd node by the boundary side of the element
587 1560 : pd_nodes[3] = fe_elem_pd_node_map.at(
588 : old_elem
589 1560 : ->neighbor_ptr(old_elem->opposite_side(
590 1560 : old_boundary_info.side_with_boundary_id(old_elem, bidit)))
591 1560 : ->id());
592 1560 : p3 = *new_mesh->node_ptr(pd_nodes[3]);
593 :
594 : // choose two neighbors of current node (total of three pd nodes) ordered in a way
595 : // such that the normal points to the fourth pd node
596 7080 : for (unsigned int i = 0; i < boundary_elem_ids.size(); ++i)
597 : {
598 5520 : Elem * nb_i = old_mesh->elem_ptr(boundary_elem_ids[i]);
599 5520 : p1 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(boundary_elem_ids[i]));
600 :
601 12720 : for (unsigned int j = i + 1; j < boundary_elem_ids.size(); ++j)
602 : {
603 7200 : Elem * nb_j = old_mesh->elem_ptr(boundary_elem_ids[j]);
604 7200 : p2 = *new_mesh->node_ptr(fe_elem_pd_node_map.at(boundary_elem_ids[j]));
605 :
606 14400 : if (old_elem->which_neighbor_am_i(nb_i) !=
607 7200 : old_elem->opposite_side(old_elem->which_neighbor_am_i(nb_j)))
608 : {
609 : should_add = true;
610 : // check whether this new elem overlaps with already created elems by the
611 : // saved edge and nodes of previously created phantom elems
612 4800 : std::pair<dof_id_type, dof_id_type> nodes_pair_i;
613 4800 : nodes_pair_i.first = std::min(eidit, boundary_elem_ids[i]);
614 4800 : nodes_pair_i.second = std::max(eidit, boundary_elem_ids[i]);
615 : if (elem_edge_nodes.count(nodes_pair_i))
616 : {
617 6375 : for (const auto & nbid : elem_edge_nodes[nodes_pair_i])
618 3425 : if (nbid != old_elem->id() && old_mesh->elem_ptr(nbid)->has_neighbor(nb_j))
619 : should_add = false;
620 : }
621 :
622 4800 : std::pair<dof_id_type, dof_id_type> nodes_pair_j;
623 4800 : nodes_pair_j.first = std::min(eidit, boundary_elem_ids[j]);
624 4800 : nodes_pair_j.second = std::max(eidit, boundary_elem_ids[j]);
625 : if (elem_edge_nodes.count(nodes_pair_j))
626 : {
627 6830 : for (const auto & nbid : elem_edge_nodes[nodes_pair_j])
628 3715 : if (nbid != old_elem->id() && old_mesh->elem_ptr(nbid)->has_neighbor(nb_i))
629 : should_add = false;
630 : }
631 :
632 4800 : if (should_add)
633 : {
634 : // check whether the normal of face formed by p0, p1 and p2 points to p3
635 : Real val =
636 2460 : ((p1(1) - p0(1)) * (p2(2) - p0(2)) - (p1(2) - p0(2)) * (p2(1) - p0(1))) *
637 2460 : (p3(0) - p0(0)) +
638 2460 : ((p1(2) - p0(2)) * (p2(0) - p0(0)) - (p1(0) - p0(0)) * (p2(2) - p0(2))) *
639 2460 : (p3(1) - p0(1)) +
640 2460 : ((p1(0) - p0(0)) * (p2(1) - p0(1)) - (p1(1) - p0(1)) * (p2(0) - p0(0))) *
641 2460 : (p3(2) - p0(2));
642 2460 : if (val > 0) // normal point to p3
643 : {
644 880 : pd_nodes[1] = fe_elem_pd_node_map.at(boundary_elem_ids[i]);
645 880 : pd_nodes[2] = fe_elem_pd_node_map.at(boundary_elem_ids[j]);
646 : }
647 : else
648 : {
649 1580 : pd_nodes[1] = fe_elem_pd_node_map.at(boundary_elem_ids[j]);
650 1580 : pd_nodes[2] = fe_elem_pd_node_map.at(boundary_elem_ids[i]);
651 : }
652 :
653 : // construct the new phantom element
654 2460 : Elem * new_elem = new Tet4;
655 2460 : new_elem->set_id(new_elem_id);
656 :
657 2460 : if (_merge_pd_blks)
658 0 : new_elem->subdomain_id() =
659 0 : min_converted_fe_blk_id + _phantom_blk_offset_number;
660 : else
661 2460 : new_elem->subdomain_id() =
662 2460 : old_elem->subdomain_id() + _phantom_blk_offset_number;
663 :
664 2460 : new_elem = new_mesh->add_elem(new_elem);
665 2460 : new_elem->set_node(0, new_mesh->node_ptr(pd_nodes[0]));
666 2460 : new_elem->set_node(1, new_mesh->node_ptr(pd_nodes[1]));
667 2460 : new_elem->set_node(2, new_mesh->node_ptr(pd_nodes[2]));
668 2460 : new_elem->set_node(3, new_mesh->node_ptr(pd_nodes[3]));
669 :
670 : // save the edges and nodes used in the new phantom elem, which will be used
671 : // for creating new phantom elements
672 2460 : elem_edge_nodes[nodes_pair_i].insert(boundary_elem_ids[j]);
673 2460 : elem_edge_nodes[nodes_pair_j].insert(boundary_elem_ids[i]);
674 :
675 2460 : ++new_elem_id;
676 :
677 2460 : new_boundary_info.add_side(new_elem, 0, _pd_blk_offset_number + bidit);
678 2460 : if (old_boundary_info.get_sideset_name(bidit) != "")
679 : {
680 0 : new_boundary_info.sideset_name(_pd_blk_offset_number + bidit) =
681 0 : "pd_side_" + old_boundary_info.get_sideset_name(bidit);
682 : }
683 : }
684 : }
685 : }
686 : }
687 1560 : }
688 : else
689 0 : mooseError("Element type ",
690 0 : old_elem->type(),
691 : " is not supported for PD sidesets construction!");
692 : }
693 1890 : }
694 : }
695 : }
696 : }
697 :
698 : // next, save non-converted or retained FE elements, if any, to new mesh
699 : std::map<dof_id_type, dof_id_type> fe_elems_map; // IDs of the same elem in the old and new meshes
700 1570 : for (const auto & eid : fe_elems)
701 : {
702 1260 : Elem * old_elem = old_mesh->elem_ptr(eid);
703 1260 : Elem * new_elem = Elem::build(old_elem->type()).release();
704 1260 : new_elem->set_id(new_elem_id);
705 1260 : new_elem->subdomain_id() = old_elem->subdomain_id();
706 1260 : new_elem = new_mesh->add_elem(new_elem);
707 5040 : for (unsigned int j = 0; j < old_elem->n_nodes(); ++j)
708 3780 : new_elem->set_node(j, new_mesh->node_ptr(fe_nodes_map.at(old_elem->node_ptr(j)->id())));
709 :
710 1260 : fe_elems_map.insert(std::make_pair(old_elem->id(), new_elem_id));
711 :
712 1260 : ++new_elem_id;
713 : }
714 :
715 : // STEP 5: convert old boundary_info to new boundary_info for retained FE mesh
716 :
717 : // peridynamics ONLY accept nodesets and sidesets
718 : // nodeset consisting single node in the converted FE block will no longer be available
719 310 : if (old_boundary_info.n_edge_conds())
720 0 : mooseError("PeridynamicsMesh does not support edgesets!");
721 :
722 : // check the existence of nodesets, if exist, build sidesets in case this hasn't been done yet
723 310 : if (old_boundary_info.n_nodeset_conds())
724 560 : old_boundary_info.build_side_list_from_node_list();
725 :
726 : // first, create a tuple to collect all sidesets (including those converted from nodesets) in the
727 : // old mesh
728 310 : auto old_fe_sbc_tuples = old_boundary_info.build_side_list();
729 : // 0: element ID, 1: side ID, 2: boundary ID
730 : // map of set of elem IDs connected to each boundary in the old mesh
731 : std::map<boundary_id_type, std::set<dof_id_type>> old_fe_bnd_elem_ids;
732 : // map of set of side ID for each elem in the old mesh
733 : std::map<dof_id_type, std::map<boundary_id_type, dof_id_type>> old_fe_elem_bnd_side_ids;
734 10435 : for (const auto & sbct : old_fe_sbc_tuples)
735 : {
736 10125 : old_fe_bnd_elem_ids[std::get<2>(sbct)].insert(std::get<0>(sbct));
737 10125 : old_fe_elem_bnd_side_ids[std::get<0>(sbct)].insert(
738 10125 : std::make_pair(std::get<2>(sbct), std::get<1>(sbct)));
739 : }
740 :
741 : // next, convert element lists in old mesh to PD nodesets in new mesh
742 : std::set<boundary_id_type> old_side_bid(old_boundary_info.get_side_boundary_ids());
743 :
744 : // loop through all old FE _sideset_ boundaries
745 1475 : for (const auto & sbid : old_side_bid)
746 11290 : for (const auto & beid : old_fe_bnd_elem_ids[sbid])
747 : if (elems_to_pd.count(beid)) // for converted FE mesh
748 : {
749 : // save corresponding boundaries on converted FE mesh to PD nodes
750 9925 : new_boundary_info.add_node(new_mesh->node_ptr(fe_elem_pd_node_map.at(beid)),
751 : sbid + _pd_blk_offset_number);
752 9925 : if (old_boundary_info.get_sideset_name(sbid) != "")
753 7280 : new_boundary_info.nodeset_name(sbid + _pd_blk_offset_number) =
754 21840 : "pd_nodes_" + old_boundary_info.get_sideset_name(sbid);
755 :
756 9925 : if (_retain_fe_mesh) // if retained, copy the corresponding boundaries, if any, to new mesh
757 : // from old mesh
758 : {
759 600 : new_boundary_info.add_side(
760 300 : new_mesh->elem_ptr(fe_elems_map[beid]), old_fe_elem_bnd_side_ids[beid][sbid], sbid);
761 300 : new_boundary_info.sideset_name(sbid) = old_boundary_info.get_sideset_name(sbid);
762 : }
763 : }
764 : else // for non-converted FE mesh, if any, copy the corresponding boundaries to new mesh
765 : // from old mesh
766 : {
767 400 : new_boundary_info.add_side(
768 200 : new_mesh->elem_ptr(fe_elems_map[beid]), old_fe_elem_bnd_side_ids[beid][sbid], sbid);
769 200 : new_boundary_info.sideset_name(sbid) = old_boundary_info.get_sideset_name(sbid);
770 : }
771 :
772 : // similar for sideset above, save _nodesets_ of non-converted and/or retained FE mesh, if any, to
773 : // new mesh
774 310 : auto old_node_bc_tuples = old_boundary_info.build_node_list();
775 : // 0: node ID, 1: boundary ID
776 : std::map<boundary_id_type, std::set<dof_id_type>> old_bnd_node_ids;
777 9800 : for (const auto & nbct : old_node_bc_tuples)
778 9490 : old_bnd_node_ids[std::get<1>(nbct)].insert(std::get<0>(nbct));
779 :
780 : std::set<boundary_id_type> old_node_bid(old_boundary_info.get_node_boundary_ids());
781 :
782 1410 : for (const auto & nbid : old_node_bid)
783 10590 : for (const auto & bnid : old_bnd_node_ids[nbid])
784 : if (fe_nodes.count(bnid))
785 : {
786 500 : new_boundary_info.add_node(new_mesh->node_ptr(fe_nodes_map.at(bnid)), nbid);
787 500 : new_boundary_info.nodeset_name(nbid) = old_boundary_info.get_sideset_name(nbid);
788 : }
789 :
790 : // create nodesets to include all PD nodes for PD blocks in the new mesh
791 49110 : for (unsigned int i = 0; i < n_pd_nodes; ++i)
792 : {
793 48800 : if (_blks_to_pd.size() > 1 && !_merge_pd_blks)
794 : {
795 : unsigned int j = 0;
796 2130 : for (const auto & blk_id : _blks_to_pd)
797 : {
798 1420 : ++j;
799 1420 : SubdomainID real_blk_id =
800 1420 : blk_id + _pd_blk_offset_number; // account for the _pd_blk_offset_number increment after
801 : // converting to PD mesh
802 1420 : if (pd_mesh.getNodeBlockID(i) == real_blk_id)
803 : {
804 710 : new_boundary_info.add_node(new_mesh->node_ptr(i), _pd_nodeset_offset_number - j);
805 710 : new_boundary_info.nodeset_name(_pd_nodeset_offset_number - j) =
806 2130 : "pd_nodes_block_" + Moose::stringify(blk_id + _pd_blk_offset_number);
807 : }
808 : }
809 : }
810 :
811 48800 : new_boundary_info.add_node(new_mesh->node_ptr(i), _pd_nodeset_offset_number);
812 48800 : new_boundary_info.nodeset_name(_pd_nodeset_offset_number) = "pd_nodes_all";
813 : }
814 :
815 : old_mesh.reset(); // destroy the old_mesh unique_ptr
816 :
817 620 : return dynamic_pointer_cast<MeshBase>(new_mesh);
818 1240 : }
|
