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 "PeridynamicsMesh.h"
11 :
12 : #include "libmesh/elem.h"
13 :
14 : registerMooseObject("PeridynamicsApp", PeridynamicsMesh);
15 :
16 : InputParameters
17 882 : PeridynamicsMesh::validParams()
18 : {
19 882 : InputParameters params = MooseMesh::validParams();
20 882 : params.addClassDescription("Mesh class to store and return peridynamics specific mesh data");
21 :
22 1764 : params.addParam<Real>("horizon_radius", "Value of horizon size in terms of radius");
23 1764 : params.addParam<Real>("horizon_number",
24 : "The material points spacing number, i.e. ratio of horizon radius to the "
25 : "effective average spacing");
26 1764 : params.addParam<Real>("bond_associated_horizon_ratio",
27 1764 : 1.5,
28 : "Ratio of bond-associated horizon to nodal horizon. This is the only "
29 : "parameters to control the size of bond-associated horizon");
30 1764 : params.addParam<std::vector<Point>>("cracks_start",
31 : "Cartesian coordinates where predefined line cracks start");
32 1764 : params.addParam<std::vector<Point>>("cracks_end",
33 : "Cartesian coordinates where predefined line cracks end");
34 1764 : params.addParam<std::vector<Real>>("cracks_width", "Widths of predefined line cracks");
35 :
36 882 : params.set<bool>("_mesh_generator_mesh") = true;
37 :
38 882 : return params;
39 0 : }
40 :
41 317 : PeridynamicsMesh::PeridynamicsMesh(const InputParameters & parameters)
42 : : MooseMesh(parameters),
43 317 : _horizon_radius(isParamValid("horizon_radius") ? getParam<Real>("horizon_radius") : 0),
44 634 : _has_horizon_number(isParamValid("horizon_number")),
45 634 : _horizon_number(_has_horizon_number ? getParam<Real>("horizon_number") : 0),
46 634 : _bah_ratio(getParam<Real>("bond_associated_horizon_ratio")),
47 1208 : _has_cracks(isParamValid("cracks_start") || isParamValid("cracks_end")),
48 634 : _dim(declareRestartableData<unsigned int>("dim")),
49 634 : _n_pdnodes(declareRestartableData<unsigned int>("n_pdnodes")),
50 634 : _n_pdbonds(declareRestartableData<unsigned int>("n_pdbonds")),
51 634 : _pdnode_average_spacing(declareRestartableData<std::vector<Real>>("pdnode_average_spacing")),
52 634 : _pdnode_horizon_radius(declareRestartableData<std::vector<Real>>("pdnode_horizon_radius")),
53 634 : _pdnode_vol(declareRestartableData<std::vector<Real>>("pdnode_vol")),
54 634 : _pdnode_horizon_vol(declareRestartableData<std::vector<Real>>("pdnode_horizon_vol")),
55 634 : _pdnode_blockID(declareRestartableData<std::vector<SubdomainID>>("pdnode_blockID")),
56 634 : _pdnode_elemID(declareRestartableData<std::vector<dof_id_type>>("pdnode_elemID")),
57 317 : _pdnode_neighbors(
58 317 : declareRestartableData<std::vector<std::vector<dof_id_type>>>("pdnode_neighbors")),
59 634 : _pdnode_bonds(declareRestartableData<std::vector<std::vector<dof_id_type>>>("pdnode_bonds")),
60 317 : _dg_neighbors(
61 317 : declareRestartableData<std::vector<std::vector<std::vector<dof_id_type>>>>("dg_neighbors")),
62 634 : _pdnode_sub_vol(declareRestartableData<std::vector<std::vector<Real>>>("pdnode_sub_vol")),
63 634 : _pdnode_sub_vol_sum(declareRestartableData<std::vector<Real>>("pdnode_sub_vol_sum")),
64 317 : _boundary_node_offset(
65 317 : declareRestartableData<std::map<dof_id_type, Real>>("boundary_node_offset")),
66 951 : _pdnode_weight_normalizer(declareRestartableData<std::vector<Real>>("pdnode_weight_normalizer"))
67 : {
68 951 : if (!(isParamValid("horizon_radius") || _has_horizon_number))
69 0 : mooseError("Must specify either horizon_radius or horizon_number to determine horizon size in "
70 : "the mesh block!");
71 :
72 317 : if (_has_cracks)
73 : {
74 60 : _cracks_start = getParam<std::vector<Point>>("cracks_start");
75 90 : _cracks_end = getParam<std::vector<Point>>("cracks_end");
76 : }
77 : else
78 : {
79 287 : _cracks_start.push_back(Point(0., 0., 0.));
80 287 : _cracks_end.push_back(Point(0., 0., 0.));
81 : }
82 :
83 317 : if (_cracks_start.size() != _cracks_end.size())
84 0 : mooseError(
85 : "Number of cracks starting points is NOT the same as number of cracks ending points!");
86 :
87 317 : if (_has_cracks)
88 : {
89 60 : if (isParamValid("cracks_width"))
90 : {
91 0 : _cracks_width = getParam<std::vector<Real>>("cracks_width");
92 0 : if (_cracks_width.size() != _cracks_start.size())
93 0 : mooseError("Number of cracks width is NOT the same as number of cracks!");
94 : }
95 : else
96 : {
97 65 : for (unsigned int i = 0; i < _cracks_start.size(); ++i)
98 35 : _cracks_width.push_back(0);
99 : }
100 : }
101 : else
102 287 : _cracks_width.push_back(0);
103 317 : }
104 :
105 : std::unique_ptr<MooseMesh>
106 248 : PeridynamicsMesh::safeClone() const
107 : {
108 248 : return _app.getFactory().copyConstruct(*this);
109 : }
110 :
111 : void
112 311 : PeridynamicsMesh::buildMesh()
113 : {
114 311 : if (!hasMeshBase())
115 : {
116 0 : auto & entry = _app.getMeshGeneratorSystem();
117 0 : _mesh = entry.getSavedMesh(entry.mainMeshGeneratorName());
118 : }
119 : _mesh->allow_renumbering(false);
120 311 : _mesh->prepare_for_use();
121 : _mesh->allow_renumbering(true);
122 311 : }
123 :
124 : unsigned int
125 182617 : PeridynamicsMesh::dimension() const
126 : {
127 182617 : return _dim;
128 : }
129 :
130 : dof_id_type
131 311 : PeridynamicsMesh::nPDNodes() const
132 : {
133 311 : return _n_pdnodes;
134 : }
135 :
136 : dof_id_type
137 311 : PeridynamicsMesh::nPDBonds() const
138 : {
139 311 : return _n_pdbonds;
140 : }
141 :
142 : void
143 311 : PeridynamicsMesh::createPeridynamicsMeshData(
144 : MeshBase & fe_mesh,
145 : std::set<dof_id_type> converted_elem_id,
146 : std::multimap<SubdomainID, SubdomainID> bonding_block_pairs,
147 : std::multimap<SubdomainID, SubdomainID> non_bonding_block_pairs)
148 : {
149 311 : _dim = fe_mesh.mesh_dimension();
150 311 : _n_pdnodes = converted_elem_id.size();
151 :
152 : // initialize data size
153 311 : _pdnode_coord.resize(_n_pdnodes);
154 311 : _pdnode_average_spacing.resize(_n_pdnodes);
155 311 : _pdnode_horizon_radius.resize(_n_pdnodes);
156 311 : _pdnode_vol.resize(_n_pdnodes);
157 311 : _pdnode_horizon_vol.resize(_n_pdnodes, 0.0);
158 311 : _pdnode_blockID.resize(_n_pdnodes);
159 311 : _pdnode_elemID.resize(_n_pdnodes);
160 311 : _pdnode_neighbors.resize(_n_pdnodes);
161 311 : _pdnode_bonds.resize(_n_pdnodes);
162 311 : _dg_neighbors.resize(_n_pdnodes);
163 311 : _pdnode_sub_vol.resize(_n_pdnodes);
164 311 : _pdnode_sub_vol_sum.resize(_n_pdnodes, 0.0);
165 311 : _pdnode_weight_normalizer.resize(_n_pdnodes, 0.0);
166 :
167 : // loop through converted fe elements to generate PD nodes structure
168 : unsigned int id = 0; // make pd nodes start at 0 in the new mesh
169 49138 : for (const auto & eid : converted_elem_id)
170 : {
171 48827 : Elem * fe_elem = fe_mesh.elem_ptr(eid);
172 : // calculate the nodes spacing as average distance between fe element with its neighbors
173 : unsigned int n_fe_neighbors = 0;
174 : Real dist_sum = 0.0;
175 259519 : for (unsigned int j = 0; j < fe_elem->n_neighbors(); ++j)
176 210692 : if (fe_elem->neighbor_ptr(j) != nullptr)
177 : {
178 196648 : dist_sum += (fe_elem->vertex_average() - fe_elem->neighbor_ptr(j)->vertex_average()).norm();
179 196648 : n_fe_neighbors++;
180 : }
181 : else // this side is on boundary and calculate the distance to the centroid
182 : {
183 : Real dist = 0.0;
184 14044 : std::vector<unsigned int> nid = fe_elem->nodes_on_side(j);
185 14044 : Point p0 = fe_elem->vertex_average();
186 14044 : Point p1 = fe_elem->point(nid[0]);
187 14044 : if (fe_elem->dim() == 2) // 2D elems
188 : {
189 8800 : Point p2 = fe_elem->point(nid.back());
190 8800 : Real area = 0.5 * std::abs(p0(0) * (p1(1) - p2(1)) + p1(0) * (p2(1) - p0(1)) +
191 8800 : p2(0) * (p0(1) - p1(1)));
192 8800 : dist = 2.0 * area / fe_elem->length(nid[0], nid[1]);
193 : }
194 : else // 3D elems
195 : {
196 5244 : Point p2 = fe_elem->point(nid[1]);
197 5244 : Point p3 = fe_elem->point(nid.back());
198 : Point vec0 = p1 - p2;
199 : Point vec1 = p1 - p3;
200 5244 : Point normal = vec0.cross(vec1);
201 5244 : normal /= normal.norm();
202 5244 : dist = std::abs(normal(0) * (p0(0) - p1(0)) + normal(1) * (p0(1) - p1(1)) +
203 5244 : normal(2) * (p0(2) - p1(2)));
204 : }
205 14044 : _boundary_node_offset.insert(std::make_pair(id, -dist));
206 14044 : }
207 :
208 48827 : _pdnode_coord[id] = fe_elem->vertex_average();
209 48827 : _pdnode_average_spacing[id] = dist_sum / n_fe_neighbors;
210 48827 : _pdnode_horizon_radius[id] =
211 48827 : (_has_horizon_number ? _horizon_number * dist_sum / n_fe_neighbors : _horizon_radius);
212 : // NOTE: PeridynamicsMesh does not support RZ/RSpherical so using volume from libmesh elem is
213 : // fine
214 48827 : _pdnode_vol[id] = fe_elem->volume();
215 48827 : _pdnode_blockID[id] = fe_elem->subdomain_id() + 1000; // set new subdomain id for PD mesh in
216 : // case FE mesh is retained
217 48827 : _pdnode_elemID[id] = fe_elem->id();
218 :
219 48827 : ++id;
220 : }
221 :
222 : // search node neighbors and create other nodal data
223 622 : createNodeHorizBasedData(bonding_block_pairs, non_bonding_block_pairs);
224 :
225 311 : createNeighborHorizonBasedData(); // applies to non-ordinary state-based model only.
226 :
227 : // total number of peridynamic bonds
228 311 : _n_pdbonds = 0;
229 49138 : for (unsigned int i = 0; i < _n_pdnodes; ++i)
230 48827 : _n_pdbonds += _pdnode_neighbors[i].size();
231 311 : _n_pdbonds /= 2;
232 :
233 : unsigned int k = 0;
234 49138 : for (unsigned int i = 0; i < _n_pdnodes; ++i)
235 2094967 : for (unsigned int j = 0; j < _pdnode_neighbors[i].size(); ++j)
236 2046140 : if (_pdnode_neighbors[i][j] > i)
237 : {
238 : // build the bond list for each node
239 1023070 : _pdnode_bonds[i].push_back(k);
240 1023070 : _pdnode_bonds[_pdnode_neighbors[i][j]].push_back(k);
241 1023070 : ++k;
242 : }
243 311 : }
244 :
245 : void
246 311 : PeridynamicsMesh::createNodeHorizBasedData(
247 : std::multimap<SubdomainID, SubdomainID> bonding_block_pairs,
248 : std::multimap<SubdomainID, SubdomainID> non_bonding_block_pairs)
249 : {
250 : // search neighbors
251 49138 : for (unsigned int i = 0; i < _n_pdnodes; ++i)
252 : {
253 : Real dis = 0.0;
254 77830576 : for (unsigned int j = 0; j < _n_pdnodes; ++j)
255 : {
256 77781749 : dis = (_pdnode_coord[i] - _pdnode_coord[j]).norm();
257 77781749 : if (dis <= 1.0001 * _pdnode_horizon_radius[i] && j != i)
258 : {
259 : bool is_interface = false;
260 1834890 : if (!bonding_block_pairs.empty())
261 : is_interface =
262 143280 : checkInterface(_pdnode_blockID[i], _pdnode_blockID[j], bonding_block_pairs);
263 :
264 1834890 : if (!non_bonding_block_pairs.empty())
265 0 : is_interface =
266 0 : !checkInterface(_pdnode_blockID[i], _pdnode_blockID[j], non_bonding_block_pairs);
267 :
268 1834890 : if (_pdnode_blockID[i] == _pdnode_blockID[j] || is_interface)
269 : {
270 : // check whether pdnode i falls in the region whose bonds may need to be removed due to
271 : // the pre-existing cracks
272 : bool intersect = false;
273 3995220 : for (unsigned int k = 0; k < _cracks_start.size(); ++k)
274 : {
275 2163700 : if (checkPointInsideRectangle(_pdnode_coord[i],
276 : _cracks_start[k],
277 : _cracks_end[k],
278 : _cracks_width[k] + 4.0 * _pdnode_horizon_radius[i],
279 2163700 : 4.0 * _pdnode_horizon_radius[i]))
280 147760 : intersect = intersect || checkCrackIntersectBond(_cracks_start[k],
281 : _cracks_end[k],
282 : _cracks_width[k],
283 : _pdnode_coord[i],
284 : _pdnode_coord[j]);
285 : }
286 : // remove bonds cross the crack to form crack surface
287 1831520 : if (!intersect)
288 : {
289 : // Use the addition balance scheme to remove unbalanced interactions
290 : // check whether j was already considered as a neighbor of i, if not, add j to i's
291 : // neighborlist
292 1819660 : if (std::find(_pdnode_neighbors[i].begin(), _pdnode_neighbors[i].end(), j) ==
293 1819660 : _pdnode_neighbors[i].end())
294 : {
295 1023070 : _pdnode_neighbors[i].push_back(j);
296 1023070 : _pdnode_horizon_vol[i] += _pdnode_vol[j];
297 1023070 : _pdnode_weight_normalizer[i] += _pdnode_horizon_radius[i] / dis;
298 : }
299 : // check whether i was also considered as a neighbor of j, if not, add i to j's
300 : // neighborlist
301 1819660 : if (std::find(_pdnode_neighbors[j].begin(), _pdnode_neighbors[j].end(), i) ==
302 1819660 : _pdnode_neighbors[j].end())
303 : {
304 1023070 : _pdnode_neighbors[j].push_back(i);
305 1023070 : _pdnode_horizon_vol[j] += _pdnode_vol[i];
306 1023070 : _pdnode_weight_normalizer[j] += _pdnode_horizon_radius[j] / dis;
307 : }
308 : }
309 : }
310 : }
311 : }
312 : }
313 311 : }
314 :
315 : bool
316 71640 : PeridynamicsMesh::checkInterface(SubdomainID pdnode_blockID_i,
317 : SubdomainID pdnode_blockID_j,
318 : std::multimap<SubdomainID, SubdomainID> blockID_pairs)
319 : {
320 : bool is_interface = false;
321 : std::pair<std::multimap<SubdomainID, SubdomainID>::iterator,
322 : std::multimap<SubdomainID, SubdomainID>::iterator>
323 : ret;
324 : // check existence of the case when i is the key and j is the value
325 : ret = blockID_pairs.equal_range(pdnode_blockID_i);
326 130550 : for (std::multimap<SubdomainID, SubdomainID>::iterator it = ret.first; it != ret.second; ++it)
327 58910 : if (pdnode_blockID_j == it->second)
328 : is_interface = true;
329 :
330 : // check existence of the case when j is the key and i is the value
331 : ret = blockID_pairs.equal_range(pdnode_blockID_j);
332 132400 : for (std::multimap<SubdomainID, SubdomainID>::iterator it = ret.first; it != ret.second; ++it)
333 60760 : if (pdnode_blockID_i == it->second)
334 : is_interface = true;
335 :
336 71640 : return is_interface;
337 : }
338 :
339 : void
340 311 : PeridynamicsMesh::createNeighborHorizonBasedData()
341 : {
342 49138 : for (unsigned int i = 0; i < _n_pdnodes; ++i)
343 : {
344 48827 : std::vector<dof_id_type> n_pd_neighbors = _pdnode_neighbors[i];
345 48827 : _dg_neighbors[i].resize(n_pd_neighbors.size());
346 48827 : _pdnode_sub_vol[i].resize(n_pd_neighbors.size(), 0.0);
347 :
348 2094967 : for (unsigned int j = 0; j < n_pd_neighbors.size(); ++j)
349 98680699 : for (unsigned int k = j; k < n_pd_neighbors.size(); ++k) // only search greater number index
350 96634559 : if ((_pdnode_coord[n_pd_neighbors[j]] - _pdnode_coord[n_pd_neighbors[k]]).norm() <=
351 96634559 : _bah_ratio * _pdnode_horizon_radius[i])
352 : {
353 : // only save the corresponding index in neighbor list, rather than the actual node id
354 : // for neighbor j
355 71176771 : _dg_neighbors[i][j].push_back(k);
356 71176771 : _pdnode_sub_vol[i][j] += _pdnode_vol[n_pd_neighbors[k]];
357 71176771 : _pdnode_sub_vol_sum[i] += _pdnode_vol[n_pd_neighbors[k]];
358 : // for neighbor k
359 71176771 : if (k > j)
360 : {
361 69130631 : _dg_neighbors[i][k].push_back(j);
362 69130631 : _pdnode_sub_vol[i][k] += _pdnode_vol[n_pd_neighbors[j]];
363 69130631 : _pdnode_sub_vol_sum[i] += _pdnode_vol[n_pd_neighbors[j]];
364 : }
365 : }
366 48827 : }
367 311 : }
368 :
369 : std::vector<dof_id_type>
370 5872443 : PeridynamicsMesh::getNeighbors(dof_id_type node_id)
371 : {
372 5872443 : return _pdnode_neighbors[node_id];
373 : }
374 :
375 : dof_id_type
376 25961336 : PeridynamicsMesh::getNeighborIndex(dof_id_type node_i, dof_id_type node_j)
377 : {
378 25961336 : std::vector<dof_id_type> n_pd_neighbors = _pdnode_neighbors[node_i];
379 25961336 : auto it = std::find(n_pd_neighbors.begin(), n_pd_neighbors.end(), node_j);
380 25961336 : if (it != n_pd_neighbors.end())
381 25961336 : return it - n_pd_neighbors.begin();
382 : else
383 0 : mooseError(
384 : "Material point ", node_j, " is not in the neighbor list of material point ", node_i);
385 :
386 : return -1;
387 25961336 : }
388 :
389 : std::vector<dof_id_type>
390 8998310 : PeridynamicsMesh::getBonds(dof_id_type node_id)
391 : {
392 8998310 : if (node_id > _n_pdnodes)
393 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
394 :
395 8998310 : return _pdnode_bonds[node_id];
396 : }
397 :
398 : std::vector<dof_id_type>
399 3200238 : PeridynamicsMesh::getBondDeformationGradientNeighbors(dof_id_type node_id, dof_id_type neighbor_id)
400 : {
401 3200238 : if (node_id > _n_pdnodes)
402 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
403 :
404 3200238 : if (neighbor_id > _pdnode_neighbors[node_id].size() - 1)
405 0 : mooseError("Querying neighbor index exceeds the available neighbors!");
406 :
407 3200238 : std::vector<dof_id_type> dg_neighbors = _dg_neighbors[node_id][neighbor_id];
408 3200238 : if (dg_neighbors.size() < _dim)
409 0 : mooseError("Not enough number of neighbors to calculate deformation gradient at PD node: ",
410 : node_id);
411 :
412 3200238 : return dg_neighbors;
413 0 : }
414 :
415 : SubdomainID
416 2047560 : PeridynamicsMesh::getNodeBlockID(dof_id_type node_id)
417 : {
418 2047560 : if (node_id > _n_pdnodes)
419 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
420 :
421 2047560 : return _pdnode_blockID[node_id];
422 : }
423 :
424 : void
425 5 : PeridynamicsMesh::setNodeBlockID(SubdomainID id)
426 : {
427 5 : _pdnode_blockID.assign(_n_pdnodes, id);
428 5 : }
429 :
430 : Point
431 202409222 : PeridynamicsMesh::getNodeCoord(dof_id_type node_id)
432 : {
433 202409222 : if (node_id > _n_pdnodes)
434 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
435 :
436 202409222 : return _pdnode_coord[node_id];
437 : }
438 :
439 : std::vector<dof_id_type>
440 0 : PeridynamicsMesh::getPDNodeIDToFEElemIDMap()
441 : {
442 0 : return _pdnode_elemID;
443 : }
444 :
445 : Real
446 106142691 : PeridynamicsMesh::getNodeVolume(dof_id_type node_id)
447 : {
448 106142691 : if (node_id > _n_pdnodes)
449 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
450 :
451 106142691 : return _pdnode_vol[node_id];
452 : }
453 :
454 : Real
455 20111786 : PeridynamicsMesh::getHorizonVolume(dof_id_type node_id)
456 : {
457 20111786 : if (node_id > _n_pdnodes)
458 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
459 :
460 20111786 : return _pdnode_horizon_vol[node_id];
461 : }
462 :
463 : Real
464 2732580 : PeridynamicsMesh::getHorizonSubsetVolume(dof_id_type node_id, dof_id_type neighbor_id)
465 : {
466 2732580 : if (node_id > _n_pdnodes)
467 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
468 :
469 2732580 : if (neighbor_id > _pdnode_neighbors[node_id].size() - 1)
470 0 : mooseError("Querying neighbor index exceeds the available neighbors!");
471 :
472 2732580 : return _pdnode_sub_vol[node_id][neighbor_id];
473 : }
474 :
475 : Real
476 2732580 : PeridynamicsMesh::getHorizonSubsetVolumeSum(dof_id_type node_id)
477 : {
478 2732580 : if (node_id > _n_pdnodes)
479 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
480 :
481 2732580 : return _pdnode_sub_vol_sum[node_id];
482 : }
483 :
484 : Real
485 13832444 : PeridynamicsMesh::getHorizonSubsetVolumeFraction(dof_id_type node_id, dof_id_type neighbor_id)
486 : {
487 13832444 : if (node_id > _n_pdnodes)
488 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
489 :
490 13832444 : if (neighbor_id > _pdnode_neighbors[node_id].size() - 1)
491 0 : mooseError("Querying neighbor index exceeds the available neighbors!");
492 :
493 13832444 : return _pdnode_sub_vol[node_id][neighbor_id] / _pdnode_sub_vol_sum[node_id];
494 : }
495 :
496 : Real
497 101520 : PeridynamicsMesh::getNodeAverageSpacing(dof_id_type node_id)
498 : {
499 101520 : if (node_id > _n_pdnodes)
500 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
501 :
502 101520 : return _pdnode_average_spacing[node_id];
503 : }
504 :
505 : Real
506 20219070 : PeridynamicsMesh::getHorizon(dof_id_type node_id)
507 : {
508 20219070 : if (node_id > _n_pdnodes)
509 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
510 :
511 20219070 : return _pdnode_horizon_radius[node_id];
512 : }
513 :
514 : Real
515 900 : PeridynamicsMesh::getBoundaryOffset(dof_id_type node_id)
516 : {
517 900 : if (_boundary_node_offset.count(node_id))
518 536 : return _boundary_node_offset[node_id];
519 : else
520 364 : return 0.0;
521 : }
522 :
523 : Real
524 12128892 : PeridynamicsMesh::getNeighborWeight(dof_id_type node_id, dof_id_type neighbor_id)
525 : {
526 12128892 : if (node_id > _n_pdnodes)
527 0 : mooseError("Querying node ID exceeds the available PD node IDs!");
528 :
529 12128892 : if (neighbor_id > _pdnode_neighbors[node_id].size() - 1)
530 0 : mooseError("Querying neighbor index exceeds the available neighbors!");
531 :
532 : Real dis =
533 12128892 : (_pdnode_coord[node_id] - _pdnode_coord[_pdnode_neighbors[node_id][neighbor_id]]).norm();
534 12128892 : Real result = _pdnode_horizon_radius[node_id] / dis / _pdnode_weight_normalizer[node_id];
535 :
536 12128892 : return result;
537 : }
538 :
539 : bool
540 2163700 : PeridynamicsMesh::checkPointInsideRectangle(
541 : Point point, Point rec_p1, Point rec_p2, Real rec_height, Real tol)
542 : {
543 2163700 : Real crack_length = (rec_p2 - rec_p1).norm();
544 : bool inside = crack_length;
545 :
546 2163700 : if (inside)
547 : {
548 : Real a = rec_p2(1) - rec_p1(1);
549 1021580 : Real b = rec_p1(0) - rec_p2(0);
550 1021580 : Real c = rec_p2(0) * rec_p1(1) - rec_p2(1) * rec_p1(0);
551 : inside =
552 1021580 : inside && (std::abs(a * point(0) + b * point(1) + c) / crack_length < rec_height / 2.0);
553 :
554 : a = rec_p2(0) - rec_p1(0);
555 : b = rec_p2(1) - rec_p1(1);
556 1021580 : c = 0.5 * (rec_p1(1) * rec_p1(1) - rec_p2(1) * rec_p2(1) - rec_p2(0) * rec_p2(0) +
557 1021580 : rec_p1(0) * rec_p1(0));
558 1021580 : inside = inside && (std::abs(a * point(0) + b * point(1) + c) / crack_length <
559 270560 : (tol + crack_length) / 2.0);
560 : }
561 :
562 2163700 : return inside;
563 : }
564 :
565 : bool
566 147760 : PeridynamicsMesh::checkCrackIntersectBond(
567 : Point crack_p1, Point crack_p2, Real crack_width, Point bond_p1, Point bond_p2)
568 : {
569 : bool intersect0 = false;
570 : bool intersect1 = false;
571 : bool intersect2 = false;
572 147760 : if ((crack_p2 - crack_p1).norm())
573 : {
574 147760 : intersect0 = checkSegmentIntersectSegment(crack_p1, crack_p2, bond_p1, bond_p2);
575 147760 : if (crack_width != 0.)
576 : {
577 0 : Real crack_len = (crack_p1 - crack_p2).norm();
578 0 : Real cos_crack = (crack_p2(0) - crack_p1(0)) / crack_len;
579 0 : Real sin_crack = (crack_p2(1) - crack_p1(1)) / crack_len;
580 0 : Real new_crack_p1x = crack_p1(0) - crack_width / 2.0 * sin_crack;
581 0 : Real new_crack_p1y = crack_p1(1) + crack_width / 2.0 * cos_crack;
582 0 : Real new_crack_p2x = crack_p2(0) - crack_width / 2.0 * sin_crack;
583 0 : Real new_crack_p2y = crack_p2(1) + crack_width / 2.0 * cos_crack;
584 : Point new_crack_p1 = Point(new_crack_p1x, new_crack_p1y, 0.);
585 : Point new_crack_p2 = Point(new_crack_p2x, new_crack_p2y, 0.);
586 0 : intersect1 = checkSegmentIntersectSegment(new_crack_p1, new_crack_p2, bond_p1, bond_p2);
587 0 : new_crack_p1x = crack_p1(0) + crack_width / 2.0 * sin_crack;
588 0 : new_crack_p1y = crack_p1(1) - crack_width / 2.0 * cos_crack;
589 0 : new_crack_p2x = crack_p2(0) + crack_width / 2.0 * sin_crack;
590 0 : new_crack_p2y = crack_p2(1) - crack_width / 2.0 * cos_crack;
591 0 : new_crack_p1 = Point(new_crack_p1x, new_crack_p1y, 0.);
592 0 : new_crack_p2 = Point(new_crack_p2x, new_crack_p2y, 0.);
593 0 : intersect2 = checkSegmentIntersectSegment(new_crack_p1, new_crack_p2, bond_p1, bond_p2);
594 : }
595 : }
596 :
597 147760 : return intersect0 || intersect1 || intersect2;
598 : }
599 :
600 : bool
601 147760 : PeridynamicsMesh::checkSegmentIntersectSegment(Point seg1_p1,
602 : Point seg1_p2,
603 : Point seg2_p1,
604 : Point seg2_p2)
605 : {
606 : // Fail if the segments share an end-point
607 1680 : if ((seg1_p1(0) == seg2_p1(0) && seg1_p1(1) == seg2_p1(1)) ||
608 147760 : (seg1_p2(0) == seg2_p1(0) && seg1_p2(1) == seg2_p1(1)) ||
609 295520 : (seg1_p1(0) == seg2_p2(0) && seg1_p1(1) == seg2_p2(1)) ||
610 5040 : (seg1_p2(0) == seg2_p2(0) && seg1_p2(1) == seg2_p2(1)))
611 : {
612 : return false;
613 : }
614 :
615 : // Fail if the segments intersect at a given end-point but not normal to the crack
616 147760 : if ((seg1_p1(1) - seg1_p2(1)) / (seg1_p1(0) - seg1_p2(0)) ==
617 295520 : (seg1_p1(1) - seg2_p1(1)) / (seg1_p1(0) - seg2_p1(0)) ||
618 : (seg1_p1(1) - seg1_p2(1)) / (seg1_p1(0) - seg1_p2(0)) ==
619 147760 : (seg1_p1(1) - seg2_p2(1)) / (seg1_p1(0) - seg2_p2(0)) ||
620 147760 : (seg2_p1(1) - seg2_p2(1)) / (seg2_p1(0) - seg2_p2(0)) ==
621 295520 : (seg2_p1(1) - seg1_p1(1)) / (seg2_p1(0) - seg1_p1(0)) ||
622 : (seg2_p1(1) - seg2_p2(1)) / (seg2_p1(0) - seg2_p2(0)) ==
623 146650 : (seg2_p1(1) - seg1_p2(1)) / (seg2_p1(0) - seg1_p2(0)))
624 : {
625 : Real COSseg1_seg2 = (seg1_p2 - seg1_p1) * (seg2_p2 - seg2_p1) /
626 2560 : ((seg1_p2 - seg1_p1).norm() * (seg2_p2 - seg2_p1).norm());
627 2560 : if (COSseg1_seg2 > -0.08715574 && COSseg1_seg2 < 0.08715574)
628 : return false;
629 : }
630 :
631 : // Translate the system so that point seg1_p1 is on the origin
632 147040 : seg1_p2(0) -= seg1_p1(0);
633 147040 : seg1_p2(1) -= seg1_p1(1);
634 147040 : seg2_p1(0) -= seg1_p1(0);
635 147040 : seg2_p1(1) -= seg1_p1(1);
636 147040 : seg2_p2(0) -= seg1_p1(0);
637 147040 : seg2_p2(1) -= seg1_p1(1);
638 :
639 : // Length of segment seg1_p1-seg1_p2
640 147040 : Real seg1_len = seg1_p2.norm();
641 :
642 : // Rotate the system so that point seg1_p2 is on the positive X axis
643 147040 : Real cos_seg1 = seg1_p2(0) / seg1_len;
644 147040 : Real sin_seg1 = seg1_p2(1) / seg1_len;
645 147040 : Real newX = seg2_p1(0) * cos_seg1 + seg2_p1(1) * sin_seg1;
646 147040 : seg2_p1(1) = seg2_p1(1) * cos_seg1 - seg2_p1(0) * sin_seg1;
647 : seg2_p1(0) = newX;
648 147040 : newX = seg2_p2(0) * cos_seg1 + seg2_p2(1) * sin_seg1;
649 147040 : seg2_p2(1) = seg2_p2(1) * cos_seg1 - seg2_p2(0) * sin_seg1;
650 : seg2_p2(0) = newX;
651 :
652 : // Fail if segment seg2_p1-seg2_p2 doesn't cross segment seg1_p1-seg1_p2
653 147040 : if ((seg2_p1(1) < 0. && seg2_p2(1) < 0.) || (seg2_p1(1) >= 0. && seg2_p2(1) >= 0.))
654 : return false;
655 :
656 : // Fail if segment seg2_p1-seg2_p2 crosses segment seg1_p1-seg1_p2 outside of segment
657 : // seg1_p1-seg1_p2
658 18195 : Real seg1_pos = seg2_p2(0) + (seg2_p1(0) - seg2_p2(0)) * seg2_p2(1) / (seg2_p2(1) - seg2_p1(1));
659 18195 : if (seg1_pos < 0. || seg1_pos > seg1_len)
660 : return false;
661 :
662 : return true;
663 : }
|
