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