Mechanics of Defect Evolution Library
Loading...
Searching...
No Matches
DislocationJunctionFormation.cpp
Go to the documentation of this file.
1/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
2 *
3 *
4 * MoDELib is distributed without any warranty under the
5 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
6 */
7
8#ifndef model_DislocationJunctionFormation_cpp_
9#define model_DislocationJunctionFormation_cpp_
10
11#include <DislocationDynamics/DislocationDynamicsBase.h>
12#include <utility> // for std::pair
13#include <algorithm>
14#include <vector>
15#include <Eigen/Dense>
16#include <SegmentSegmentDistance.h>
17#include <SweepPlane.h>
18#include <TypeTraits.h>
19#include <DislocationDynamics/StressStraight.h>
20
21#include <DislocationJunctionFormation.h>
22
23#ifndef NDEBUG
24#define VerboseJunctions(N,x) if(verboseJunctions>=N){std::cout<<x;}
25#else
26#define VerboseJunctions(N,x)
27#endif
28
29
30#include <DislocationNetwork.h>
31
32namespace model
33{
34
35 template <int dim>
36 void DislocationJunctionFormation<dim>::insertIntersection(std::deque<IntersectionTypeContainerType> &intersectionContainer,
37 const NetworkLinkType *const linkA,
38 const NetworkLinkType *const linkB,
39 const SegmentSegmentDistance<dim> &ssd,
40 const double &currentcCollisionTOL,
41 const bool &bndJunction, const bool &gbndJunction,const int event)
42 {
43 VerboseJunctions(2, "insertIntersection " << linkA->tag() << " [ " << linkA->source->get_P().transpose() << ", " << linkA->sink->get_P().transpose()
44 << " ]"
45 << "," << linkB->tag() << " [ " << linkB->source->get_P().transpose() << ", " << linkB->sink->get_P().transpose()
46 << " ]" << std::endl;);
47 VerboseJunctions(2, "insertIntersection "
48 << "ssd.dMin=" << ssd.dMin << std::endl;);
49 VerboseJunctions(2, "insertIntersection "
50 << "currentcCollisionTOL=" << currentcCollisionTOL << std::endl;);
51
52 if (ssd.dMin < currentcCollisionTOL)
53 {
54 bool isValidJunction(bndJunction || gbndJunction);
55 if (!isValidJunction && !linkA->isBoundarySegment() && !linkB->isBoundarySegment() && !linkA->isGrainBoundarySegment() && !linkB->isGrainBoundarySegment())
56 { // Check force condition for internal segments
57
58 const VectorDim chordA(linkA->sink->get_P() - linkA->source->get_P());
59 const double LA(chordA.norm());
60 const VectorDim chordB(linkB->sink->get_P() - linkB->source->get_P());
61 const double LB(chordB.norm());
62
63 if (LA > FLT_EPSILON && LB > FLT_EPSILON)
64 {
65 if (ssd.dMin > FLT_EPSILON)
66 {
67 StressStraight<dim> stressA(DN.ddBase.poly,ssd.x0 - infiniteLineLength / LA * chordA,
68 ssd.x0 + infiniteLineLength / LA * chordA,
69 linkA->burgers(),DN.ddBase.dislocationDynamics().EwaldLength);
70
71 StressStraight<dim> stressB(DN.ddBase.poly,ssd.x1 - infiniteLineLength / LB * chordB,
72 ssd.x1 + infiniteLineLength / LB * chordB,
73 linkB->burgers(),DN.ddBase.dislocationDynamics().EwaldLength);
74
75 VerboseJunctions(3, "Non-intersecting pair" << std::endl;);
76 const VectorDim forceOnA = (stressB.stress(ssd.x0) * linkA->burgers()).cross(chordA);
77 const VectorDim forceOnB = (stressA.stress(ssd.x1) * linkB->burgers()).cross(chordB);
78 const VectorDim dxShift(ssd.x1 - ssd.x0);
79 if (forceOnA.dot(dxShift) > FLT_EPSILON && forceOnB.dot(dxShift) < -FLT_EPSILON)
80 {
81 VerboseJunctions(3, "attractive pair 1" << std::endl;);
82 isValidJunction = true; // for non-parallel lines this neglects the energy of rotation
83 }
84 else
85 {
86 VerboseJunctions(3, "non-attractive pair" << std::endl;);
87 }
88 }
89 else
90 {
91 VerboseJunctions(3, "Intersecting pair .. Determining via frank rule" << std::endl;);
92 isValidJunction=(linkA->burgers().dot(linkB->burgers())*linkA->chord().dot(linkB->chord())<=0.0);
93 }
94 }
95 }
96
97 if (isValidJunction)
98 {
99 #ifdef _OPENMP
100 intersectionContainer[omp_get_thread_num()].emplace_back(linkA->key,
101 linkB->key,
102 ssd,event);
103 #else
104 intersectionContainer[0].emplace_back(linkA->key,
105 linkB->key,
106 ssd,event);
107 #endif
108 }
109 }
110 else
111 {
112 VerboseJunctions(3, "dMin=" << ssd.dMin << ", collisionTol=" << currentcCollisionTOL << std::endl;);
113 }
114 }
115
116 template <int dim>
117 void DislocationJunctionFormation<dim>::findIntersections(std::deque<IntersectionTypeContainerType>& intersectionContainer,
118 const size_t& nThreads)
119
120 {
123
124 const auto t0= std::chrono::system_clock::now();
125 std::cout<<" Finding collisions "<<std::flush;
126
127 // Use SweepPlane to compute possible intersections
128 SweepPlane<NetworkLinkType,dim> swp;
129 for(const auto& link : DN.networkLinks())
130 {
131 const auto sharedLink(link.second.lock());
132
133 if( (!sharedLink->hasZeroBurgers() )
134 || sharedLink->isBoundarySegment() || sharedLink->isGrainBoundarySegment() )
135
136 {
137 const double x0=sharedLink->source->get_P()(0);
138 const double x1=sharedLink->sink->get_P()(0);
139 swp.addSegment(std::min(x0,x1)-0.5*collisionTol,
140 std::max(x0,x1)+0.5*collisionTol,*sharedLink);
141 }
142 }
143 swp.computeIntersectionPairs();
144 std::cout<<"("<<swp.potentialIntersectionPairs().size()<<" sweep-line pairs) "<<defaultColor<<std::flush;
145
146
147 std::deque<std::pair<const NetworkLinkType*,const NetworkLinkType*>> reducedIntersectionPairs;
148 for(size_t k=0;k<swp.potentialIntersectionPairs().size();++k)
149 {
150 const auto& linkA(swp.potentialIntersectionPair(k).first);
151 const auto& linkB(swp.potentialIntersectionPair(k).second);
152 const VectorDim& sourceA(linkA->source->get_P());
153 const VectorDim& sinkA(linkA->sink->get_P());
154 const VectorDim& sourceB(linkB->source->get_P());
155 const VectorDim& sinkB(linkB->sink->get_P());
156
157 // SweepPlane sweeps along x direction, so now check possible overlap in y and z direction
158 bool yDontOverlap((std::min(sourceA(1),sinkA(1))>std::max(sourceB(1),sinkB(1))+collisionTol) || (std::max(sourceA(1),sinkA(1))<std::min(sourceB(1),sinkB(1))-collisionTol));
159 bool zDontOverlap((std::min(sourceA(2),sinkA(2))>std::max(sourceB(2),sinkB(2))+collisionTol) || (std::max(sourceA(2),sinkA(2))<std::min(sourceB(2),sinkB(2))-collisionTol));
160
161 if(!yDontOverlap && !zDontOverlap)
162 {
163 reducedIntersectionPairs.emplace_back(linkA,linkB);
164 }
165 }
166
167 std::cout<<" ("<<reducedIntersectionPairs.size()<<" reduced pairs) "<<magentaColor<<" ["<<(std::chrono::duration<double>(std::chrono::system_clock::now()-t0)).count()<<" sec]"<<defaultColor<<std::endl;
168
169 const auto t1= std::chrono::system_clock::now();
170 std::cout<<" Selecting junctions ("<<nThreads<<" threads): "<<std::flush;
171
173 // #ifdef _OPENMP
174 // #pragma omp parallel for
175 // #endif
176 for(size_t k=0;k<reducedIntersectionPairs.size();++k)
177 {
178
179 const auto& linkA(reducedIntersectionPairs[k].first);
180 const auto& linkB(reducedIntersectionPairs[k].second);
181
182 findPairIntersections(intersectionContainer,linkA,linkB);
183 }
184
185 int nIntersections=0;
186 for (const auto& intersectionByThreadContainer : intersectionContainer)
187 {
188 nIntersections+=intersectionByThreadContainer.size();
189 }
190 std::cout<<nIntersections<<" physical junctions "<<std::flush;
191 std::cout<<magentaColor<<" ["<<(std::chrono::duration<double>(std::chrono::system_clock::now()-t1)).count()<<" sec]"<<defaultColor<<std::endl;
192
193 }
194
195 template <int dim>
196 void DislocationJunctionFormation<dim>::findPairIntersections(
197 std::deque<IntersectionTypeContainerType>& intersectionContainer,
198 const NetworkLinkType* linkA,const NetworkLinkType* linkB)
199 {
200 VerboseJunctions(2,"checking intersection "<<linkA->tag()<<","<<linkB->tag()<<std::endl;);
201
202
203 const bool linkAisBnd(linkA->isBoundarySegment());
204 const bool linkBisBnd(linkB->isBoundarySegment());
205 const bool linkAisGBnd(linkA->isGrainBoundarySegment());
206 const bool linkBisGBnd(linkB->isGrainBoundarySegment());
207
208
209 bool bndJunction=false;
210 bool gbndJunction=false;
211 const auto facesB=linkB->meshFaces();
212 // Boundary merging requires a common physical face, not merely a
213 // boundary label on one segment or parallel faces elsewhere.
214 if(linkA->chord().norm()>FLT_EPSILON && linkB->chord().norm()>FLT_EPSILON &&
215 linkA->chord().normalized().cross(linkB->chord().normalized()).norm()<FLT_EPSILON &&
216 linkA->glidePlaneNormal().cross(linkB->glidePlaneNormal()).norm()<FLT_EPSILON)
217 {
218 for(const auto* face : linkA->meshFaces())
219 {
220 if(facesB.count(face))
221 {
222 const bool grainBoundary=face->regionIDs.first!=face->regionIDs.second ||
223 (face->periodicFacePair.second && face->regionIDs!=face->periodicFacePair.second->regionIDs);
224 gbndJunction=gbndJunction || grainBoundary;
225 bndJunction=bndJunction || (face->isExternal() && !face->periodicFacePair.second);
226 }
227 }
228 }
229
230 VerboseJunctions(2,"links are bnd: "<<linkAisBnd<<" "<<linkBisBnd<<std::endl;);
231 VerboseJunctions(2,"links are grainBnd: "<<linkAisGBnd<<" "<<linkBisGBnd<<std::endl;);
232 VerboseJunctions(2,"bndJunction: "<<bndJunction<<std::endl;);
233 VerboseJunctions(2,"gbndJunction: "<<gbndJunction<<std::endl;);
234
235 const bool intersectionIsSourceSource(linkA->source->sID==linkB->source->sID);
236 const bool intersectionIsSourceSink( linkA->source->sID==linkB-> sink->sID);
237 const bool intersectionIsSinkSource( linkA-> sink->sID==linkB->source->sID);
238 const bool intersectionIsSinkSink( linkA-> sink->sID==linkB-> sink->sID);
239
240 double currentcCollisionTOL=collisionTol;
241
242 const auto pcPlanes(linkA->parallelAndCoincidentGlidePlanes(linkB->glidePlanes()));
243 VerboseJunctions(2,"pcPlanes.size()= "<<pcPlanes.size()<<std::endl;);
244 for(const auto& pair : pcPlanes)
245 {
246 if(pair.first==pair.second)
247 {// a common coincident plane found
248 const double cosTheta=(intersectionIsSourceSource||intersectionIsSinkSink)? linkA->chord().normalized().dot(linkB->chord().normalized()) : ((intersectionIsSourceSink ||intersectionIsSinkSource)? -linkA->chord().normalized().dot(linkB->chord().normalized()) : -1.0);
249 if(cosTheta<0.7)
250 {// links are either disconnected, or connected to a node and forming a large angle at that node. Reduce tolerance
251 currentcCollisionTOL=FLT_EPSILON;
252 }
253 break;
254 }
255 else
256 {// segments on parallel (non-coincident) planes
257 currentcCollisionTOL=FLT_EPSILON;
258 }
259 }
260
261 if(intersectionIsSourceSource)
262 {
263 // intersect sink of A with link B
264 SegmentSegmentDistance<dim> ssdA(linkA->sink->get_P(), // fake degenerete segment at sink of A
265 linkA->sink->get_P(), // fake degenerete segment at sink of A
266 linkB->source->get_P(),
267 linkB->sink->get_P(),
268 1.0);
269 insertIntersection(intersectionContainer,linkA,linkB,ssdA,currentcCollisionTOL,bndJunction,gbndJunction,1);
270
271 // intersect sink of B with link A
272 SegmentSegmentDistance<dim> ssdB(linkA->source->get_P(),
273 linkA->sink->get_P(),
274 linkB->sink->get_P(), // fake degenerete segment at sink of B
275 linkB->sink->get_P(), // fake degenerete segment at sink of B
276 1.0);
277 insertIntersection(intersectionContainer,linkA,linkB,ssdB,currentcCollisionTOL,bndJunction,gbndJunction,2);
278
279 }
280 else if(intersectionIsSourceSink)
281 {
282 // intersect sink of A with link B
283 SegmentSegmentDistance<dim> ssdA(linkA->sink->get_P(), // fake degenerete segment at sink of A
284 linkA->sink->get_P(), // fake degenerete segment at sink of A
285 linkB->source->get_P(),
286 linkB->sink->get_P(),
287 1.0);
288 insertIntersection(intersectionContainer,linkA,linkB,ssdA,currentcCollisionTOL,bndJunction,gbndJunction,1);
289
290 // intersect source of B with link A
291 SegmentSegmentDistance<dim> ssdB(linkA->source->get_P(),
292 linkA->sink->get_P(),
293 linkB->source->get_P(), // fake degenerete segment at source of B
294 linkB->source->get_P(), // fake degenerete segment at source of B
295 0.0);
296 insertIntersection(intersectionContainer,linkA,linkB,ssdB,currentcCollisionTOL,bndJunction,gbndJunction,2);
297 }
298 else if(intersectionIsSinkSource)
299 {
300
301 // intersect source of A with link B
302 SegmentSegmentDistance<dim> ssdA(linkA->source->get_P(), // fake degenerete segment at source of A
303 linkA->source->get_P(), // fake degenerete segment at source of A
304 linkB->source->get_P(),
305 linkB->sink->get_P(),
306 0.0);
307 insertIntersection(intersectionContainer,linkA,linkB,ssdA,currentcCollisionTOL,bndJunction,gbndJunction,1);
308
309
310 // intersect sink of B with link A
311 SegmentSegmentDistance<dim> ssdB(linkA->source->get_P(),
312 linkA->sink->get_P(),
313 linkB->sink->get_P(), // fake degenerete segment at sink of B
314 linkB->sink->get_P(), // fake degenerete segment at sink of B
315 1.0);
316 insertIntersection(intersectionContainer,linkA,linkB,ssdB,currentcCollisionTOL,bndJunction,gbndJunction,2);
317 }
318 else if(intersectionIsSinkSink)
319 {
320 // intersect source of A with link B
321 SegmentSegmentDistance<dim> ssdA(linkA->source->get_P(), // fake degenerete segment at source of A
322 linkA->source->get_P(), // fake degenerete segment at source of A
323 linkB->source->get_P(),
324 linkB->sink->get_P(),
325 0.0);
326 insertIntersection(intersectionContainer,linkA,linkB,ssdA,currentcCollisionTOL,bndJunction,gbndJunction,1);
327
328
329 // intersect source of B with link A
330 SegmentSegmentDistance<dim> ssdB(linkA->source->get_P(),
331 linkA->sink->get_P(),
332 linkB->source->get_P(), // fake degenerete segment at source of B
333 linkB->source->get_P(), // fake degenerete segment at source of B
334 0.0);
335 insertIntersection(intersectionContainer,linkA,linkB,ssdB,currentcCollisionTOL,bndJunction,gbndJunction,2);
336
337 }
338 else
339 {// segments don't share vertices
340
341 SegmentSegmentDistance<dim> ssd(linkA->source->get_P(),
342 linkA->sink->get_P(),
343 linkB->source->get_P(),
344 linkB->sink->get_P()
345 );
346 insertIntersection(intersectionContainer,linkA,linkB,ssd,currentcCollisionTOL,bndJunction,gbndJunction);
347 }
348 }
349
350 // Yash's version(Working)
351 template <int dim>
352 std::pair<std::shared_ptr<typename DislocationJunctionFormation<dim>::NetworkNodeType>, bool> DislocationJunctionFormation<dim>::junctionNode(const double &t,
353 const VectorDim &x,
354 const std::shared_ptr<NetworkLinkType> &L)
355 {
356 // VerboseJunctions(4, "JunctionNode for segment: " << key.first << "->" << key.second << " @ " << t << std::endl;);
357
358 auto Nclose = t < 0.5 ? L->source : L->sink;
359 auto Nfar = t >= 0.5 ? L->source : L->sink;
360 VerboseJunctions(4, "Nclose=" << Nclose->sID << " (GP Size)=> " << Nclose->glidePlanes().size() << std::endl;);
361 VerboseJunctions(4, "Nfar=" << Nfar->sID << " (GP Size)=> " << Nfar->glidePlanes().size() << std::endl;);
362 // if(t>FLT_EPSILON && t<1.0-FLT_EPSILON && !Nclose->isBoundaryNode() && !Nfar->isBoundaryNode()) //This condition is the problem It needs to be enabled
363 if (t > FLT_EPSILON && t < 1.0 - FLT_EPSILON)
364 { // intersection point is not an end node
365 if ((Nclose->get_P() - x).norm() > DN.networkRemesher.Lmin)
366 { // far enough from Nclose
367 VerboseJunctions(4, "JunctionNode case a" << std::endl;);
368
369 const VectorDim snappedPosition(L->glidePlanes().size() > 0 ? L->snapToGlidePlanes(x) : x);
370 const auto newNetNode(DN.networkNodes().create(snappedPosition, L->source->get_V() * (1.0 - t) + L->sink->get_V() * t,L->source->climbVelocityScalar*(1.0-t)+L->sink->climbVelocityScalar*t,L->source->velocityReduction() * (1.0 - t) + L->sink->velocityReduction() * t));
371 DN.expandNetworkLink(L, newNetNode);
372 return std::make_pair(newNetNode, true);
373 }
374 else
375 { // close to from Nclose
376 if (Nclose->isMovableTo(x) && !Nclose->isBoundaryNode())
377 {
378 VerboseJunctions(4, "JunctionNode case b" << std::endl;);
379 return std::make_pair(Nclose, false);
380 }
381 else
382 {
383 if ((Nfar->get_P() - x).norm() > DN.networkRemesher.Lmin)
384 {
385 VerboseJunctions(4, "JunctionNode case c" << std::endl;);
386
387 const VectorDim snappedPosition(L->glidePlanes().size() > 0 ? L->snapToGlidePlanes(x) : x);
388
389 const auto newNetNode(DN.networkNodes().create(snappedPosition, L->source->get_V() * (1.0 - t) + L->sink->get_V() * t, L->source->climbVelocityScalar*(1.0-t)+L->sink->climbVelocityScalar*t,L->source->velocityReduction() * (1.0 - t) + L->sink->velocityReduction() * t));
390 // std::cout << "x is " << x.transpose() << std::endl;
391
392 // std::cout<<"Expanding network Links c "<<L->tag()<<" [ "<<L->source->get_P().transpose()<<" -> "<<L->sink->get_P().transpose()<<" ] GlidePlane size "
393 // <<" [ "<<L->source->glidePlanes().size()<<" --> "<<L->sink->glidePlanes().size() <<" ] "
394 // <<L->glidePlanes().size()<<" loopLinks size C"<<L->loopLinks().size()<<std::endl;
395
396 DN.expandNetworkLink(L, newNetNode);
397 // std::cout<<"Expanded network Links"<<std::endl;
398
399 return std::make_pair(newNetNode, true);
400 }
401 else
402 {
403 if (Nfar->isMovableTo(x) && !Nfar->isBoundaryNode())
404 {
405 VerboseJunctions(4, "JunctionNode case d" << std::endl;);
406 return std::make_pair(Nfar, false);
407 }
408 else
409 {
410 VerboseJunctions(4, "JunctionNode case e" << std::endl;);
411 // std::cout<<"x is "<<x.transpose()<<std::endl;
412 // std::cout<<"Expanding network Links E "<<L->tag()<<" [ "<<L->source->get_P().transpose()<<" -> "<<L->sink->get_P().transpose()<<" ] GlidePlane size "
413 // <<" [ "<<L->source->glidePlanes().size()<<" --> "<<L->sink->glidePlanes().size() <<" ] "
414 // <<L->glidePlanes().size()<<" loopLinks size E"<<L->loopLinks().size()<<std::endl;
415
416 const VectorDim snappedPosition(L->glidePlanes().size() > 0 ? L->snapToGlidePlanes(x) : x);
417
418 const auto newNetNode(DN.networkNodes().create(snappedPosition, L->source->get_V() * (1.0 - t) + L->sink->get_V() * t, L->source->climbVelocityScalar*(1.0-t)+L->sink->climbVelocityScalar*t,L->source->velocityReduction() * (1.0 - t) + L->sink->velocityReduction() * t));
419 DN.expandNetworkLink(L, newNetNode);
420 return std::make_pair(newNetNode, true);
421 }
422 }
423 }
424 }
425 }
426 else
427 {// junction node is an existing node
428 bool physicalEndpoint=false;
429 for(const auto* ln : Nclose->loopNodes())
430 {
431 physicalEndpoint=physicalEndpoint || ln->isPhysicalBoundaryNode();
432 }
433 if (!Nclose->isBoundaryNode() ||
434 (physicalEndpoint && (Nclose->get_P()-x).norm()<=FLT_EPSILON))
435 {
436 VerboseJunctions(4, "JunctionNode case f" << std::endl;);
437 return std::make_pair(Nclose, false);
438 }
439 else
440 {
441 VerboseJunctions(4, "JunctionNode case g (returning nullptr)" << std::endl;);
442 return std::make_pair(nullptr, false);
443 }
444 }
445 }
446
447 //New Implemnetation with update of the boundary nodes
448 template <int dim>
449 size_t DislocationJunctionFormation<dim>::contractJunctions(const std::deque<IntersectionTypeContainerType>& intersectionContainer)
450 {
451 const auto t0= std::chrono::system_clock::now();
452 std::cout<<" : "<<std::flush;
453
454 DN.danglingBoundaryLoopNodes.clear();
455 size_t nContracted=0;
456 size_t nDeferred=0;
457 std::set<size_t> changedNodes;
458
459 for (const auto& intersectionByThreadContainer : intersectionContainer)
460 {
461 for (const auto& intersection : intersectionByThreadContainer)
462 {
463 const KeyType& key1(std::get<0>(intersection));
464 const KeyType& key2(std::get<1>(intersection));
465
466 const auto L1(DN.networkLinks().get(key1));
467 const auto L2(DN.networkLinks().get(key2));
468
469 if(L1 && L2) // Links exist
470 {
471 // Inspect whole incident contours before querying derived boundary
472 // topology. A shared node makes two otherwise distinct loops
473 // dependent. IDs also survive loop splitting during contraction.
474 std::set<size_t> footprint;
475 std::set<size_t> visitedLoops;
476 for(const auto& node : {L1->source,L1->sink,L2->source,L2->sink})
477 {
478 for(const auto* occurrence : node->loopNodes())
479 {
480 if(visitedLoops.insert(occurrence->loop()->sID).second)
481 {
482 for(const auto* link : occurrence->loop()->linkSequence())
483 {
484 footprint.insert(link->source->networkNode->sID);
485 }
486 }
487 }
488 }
489 const bool dependent=std::any_of(footprint.begin(),footprint.end(),[&](const size_t id)
490 {
491 return changedNodes.count(id)!=0;
492 });
493 if(dependent)
494 {
495 ++nDeferred;
496 continue;
497 }
498 // Recheck only queued pairs, not the global sweep. Earlier
499 // contractions can move endpoints or change link membership
500 // without changing the two endpoint IDs.
501 std::deque<IntersectionTypeContainerType> refreshed(1);
502 findPairIntersections(refreshed,L1.get(),L2.get());
503 if(refreshed.front().empty())
504 {
505 continue;
506 }
507 // Keep the two endpoint-versus-segment events distinct.
508 // If connectivity changed, use the newly valid pair event.
509 auto refreshedEvent=refreshed.front().begin();
510 for(auto it=refreshed.front().begin();it!=refreshed.front().end();++it)
511 {
512 if(std::get<3>(*it)==std::get<3>(intersection))
513 {
514 refreshedEvent=it;
515 break;
516 }
517 }
518 const auto ssd=std::get<2>(*refreshedEvent);
519 const double t=ssd.t;
520 const double u=ssd.u;
521 VerboseJunctions(1,"forming junction "<<key1.first<<"->"<<key1.second<<", "
522 /* */ <<key2.first<<"->"<<key2.second<<", "
523 /* */ <<"dMin="<<ssd.dMin<<", "
524 /* */ <<"@ ("<<t<<","<<u<<"), "<<std::flush);
525 VerboseJunctions(4, "JunctionNode for segment: " << key1.first << "->" << key1.second << " @ " << t << std::endl;);
526 std::pair<std::shared_ptr<NetworkNodeType>,bool> Ni(junctionNode(t,ssd.x0,L1));
527 VerboseJunctions(4, "JunctionNode for segment: " << key2.first << "->" << key2.second << " @ " << u << std::endl;);
528 std::pair<std::shared_ptr<NetworkNodeType>,bool> Nj(junctionNode(u,ssd.x1,L2));
529
530 if (Ni.first && Nj.first)
531 {
532
533 VerboseJunctions(1, "contracting " << Ni.first->sID << " " << Nj.first->sID << std::endl;);
534
535 if (Ni.first->sID != Nj.first->sID)
536 {
537 const bool success = DN.nodeContractor.contract(Ni.first,Nj.first,false);
538 VerboseJunctions(1, "success= " << success << std::endl;);
539
540 nContracted += success;
541 if (!success)
542 {
543 if (Ni.second)
544 {
545 DN.removeNetworkNode(Ni.first->sID);
546 }
547 if (Nj.second)
548 {
549 DN.removeNetworkNode(Nj.first->sID);
550 }
551 }
552 else
553 {
554 changedNodes.insert(footprint.begin(),footprint.end());
555 changedNodes.insert(Ni.first->sID);
556 changedNodes.insert(Nj.first->sID);
557 }
558 }
559 }
560 else
561 {
562 //Either Ni or Nj might have been created...
563 //Get rid of that
564 if (Ni.first)
565 {
566 if (Ni.second)
567 {
568 DN.removeNetworkNode(Ni.first->sID);
569 }
570 }
571
572 if (Nj.first)
573 {
574 if (Nj.second)
575 {
576 DN.removeNetworkNode(Nj.first->sID);
577 }
578 }
579
580 }
581
582 }
583 }
584 }
585
586 if(nContracted || !DN.danglingBoundaryLoopNodes.empty())
587 {
588 DN.updateBoundaryNodes();
589 }
590 std::cout<<" ("<<nContracted<<" contracted, "<<nDeferred<<" deferred)"<<magentaColor
591 <<" ["<<(std::chrono::duration<double>(std::chrono::system_clock::now()-t0)).count()
592 <<" sec]"<<defaultColor<<std::endl;
593
594 return nContracted;
595 }
596
597 //This version also compares the network links burgers across the boundary
598 template <int dim>
599 void DislocationJunctionFormation<dim>::glissileJunctions(const double &dx)
600 {
601 const auto t0 = std::chrono::system_clock::now();
602 std::cout << " Forming Glissile Junctions: " << std::flush;
603
604 std::deque<std::tuple<std::shared_ptr<NetworkNodeType>, std::shared_ptr<NetworkNodeType>, size_t, size_t>> glissDeq;
605
606 for (const auto &linkIter : DN.networkLinks())
607 {
608 const auto link(linkIter.second.lock());
609
610 // if (link->isSessile() && link->loopLinks().size() > 1) // a junction
611 if ( link->loopLinks().size() > 1 // a junction
612 && link->chordLength() > dx
613 && link->isSessile()
614 )
615 {
616
617 const VectorDim unitChord(link->chord() / link->chordLength());
618
619 if (!link->isGrainBoundarySegment() && !link->isBoundarySegment())
620 {
621
622 VerboseJunctions(2, "glissele junction, segment " << link->tag() << std::endl;);
623
624 for (const auto &gr : link->grains())
625 {
626 for (size_t s = 0; s < gr->slipSystems().size(); ++s)
627 {
628 const auto &slipSystem(gr->slipSystems().at(s));
629 if ((slipSystem->s.cartesian() - link->burgers()).norm() < FLT_EPSILON && fabs(slipSystem->n.cartesian().normalized().dot(unitChord)) < FLT_EPSILON)
630 {
631 VerboseJunctions(3, "glissDeq, emplacing" << std::endl;);
632
633 glissDeq.emplace_back(link->source, link->sink, gr->grainID, s);
634 }
635 }
636 }
637 }
638 }
639 }
640
641 size_t formedJunctions=0;
642 for(const auto& [source,sink,grainID,slipID] : glissDeq)
643 {
644 const auto link=DN.networkLinks().get(std::make_pair(source->sID,sink->sID));
645 const auto& slip=DN.ddBase.poly.grain(grainID)->slipSystems().at(slipID);
646 // Earlier seeds may already have consumed this periodic continuation.
647 if(!link || link->hasZeroBurgers() || !link->isSessile() ||
648 link->isBoundarySegment() || link->isGrainBoundarySegment() ||
649 (link->burgers()-slip->s.cartesian()).norm()>FLT_EPSILON)
650 {
651 continue;
652 }
653
654 // Orient each contributing loop's walk from network source to sink.
655 // Only artificial periodic subdivisions are traversed: physical
656 // faces and original loop nodes terminate the continuation.
657 std::deque<const LoopNodeType*> path;
658 VectorDim originShift=VectorDim::Zero();
659 bool valid=true;
660 for(const auto* occurrence : link->loopLinks())
661 {
662 const bool forward=occurrence->source->networkNode==source;
663 const auto* a=(forward? occurrence->source : occurrence->sink).get();
664 const auto* b=(forward? occurrence->sink : occurrence->source).get();
665 std::deque<const LoopNodeType*> candidate{a,b};
666 std::set<const LoopNodeType*> visited{a,b};
667 for(const bool left : {true,false})
668 {
669 const auto* current=left? candidate.front() : candidate.back();
670 while(current->isPeriodicSubdivision())
671 {
672 const bool next=left? !forward : forward;
673 const auto* adjacent=next? current->next.first : current->prev.first;
674 if(!adjacent || !visited.insert(adjacent).second)
675 {
676 valid=false;
677 break;
678 }
679 if(left)
680 {
681 candidate.push_front(adjacent);
682 }
683 else
684 {
685 candidate.push_back(adjacent);
686 }
687 current=adjacent;
688 }
689 if(!valid)
690 {
691 break;
692 }
693 }
694 if(!valid)
695 {
696 break;
697 }
698 const VectorDim shift=a->periodicPlanePatch()->shift;
699 if(path.empty())
700 {
701 path=candidate;
702 originShift=shift;
703 }
704 else
705 {
706 // Sharing must agree through the whole continuation, not
707 // just at the original pair of network nodes. Do not choose
708 // an unrelated loop visiting a shared endpoint.
709 if(path.size()!=candidate.size())
710 {
711 valid=false;
712 break;
713 }
714 for(size_t k=0;k<path.size();++k)
715 {
716 if(path[k]->networkNode!=candidate[k]->networkNode ||
717 (path[k]->periodicPlanePatch()->shift-originShift-
718 candidate[k]->periodicPlanePatch()->shift+shift).norm()>FLT_EPSILON)
719 {
720 valid=false;
721 break;
722 }
723 }
724 }
725 if(!valid)
726 {
727 break;
728 }
729 }
730 if(!valid || path.size()<2)
731 {
732 continue;
733 }
734
735 const VectorDim midpoint=0.5*(source->get_P()+sink->get_P());
736 const auto glidePlane=DN.ddBase.dislocationDynamics().glidePlaneFactory.getFromKey(GlidePlaneKey<dim>(midpoint,slip->n));
737 const auto periodicPlane=DN.ddBase.dislocationDynamics().periodicGlidePlaneFactory.get(glidePlane->key);
738 using EdgePair=std::pair<std::shared_ptr<PeriodicPlaneEdge<dim>>,std::shared_ptr<PeriodicPlaneEdge<dim>>>;
739 struct SeedPoint
740 {
741 std::shared_ptr<NetworkNodeType> node;
742 VectorDim position;
743 std::shared_ptr<PeriodicPlanePatch<dim>> patch;
744 EdgePair edges;
745 };
746 std::vector<SeedPoint> prepared;
747 std::set<const NetworkLinkType*> physicalLinks;
748 size_t midpointAfter=path.size();
749 for(size_t k=0;k<path.size();++k)
750 {
751 const auto* node=path[k];
752 const VectorDim shift=node->periodicPlanePatch()->shift-originShift;
753 const VectorDim position=node->networkNode->get_P()-shift;
754 if(std::abs((position-glidePlane->P).dot(glidePlane->unitNormal))>FLT_EPSILON)
755 {
756 valid=false;
757 break;
758 }
759 const auto patch=periodicPlane->getPatch(shift);
760 EdgePair edges;
761 for(const auto& edge : patch->edges())
762 {
763 // Geometry assigns metadata on the target plane; it never
764 // chooses the continuation. Corners retain both patch edges.
765 if(edge->meshIntersection->contains(node->networkNode->get_P()))
766 {
767 if(!edges.first)
768 {
769 edges.first=edge;
770 }
771 else if(!edges.second)
772 {
773 edges.second=edge;
774 }
775 else
776 {
777 valid=false;
778 }
779 }
780 }
781 if(node->isPeriodicSubdivision() && !edges.first)
782 {
783 valid=false;
784 }
785 prepared.push_back({node->networkNode,position,patch,edges});
786 if(k+1<path.size())
787 {
788 const auto* next=path[k+1];
789 const auto* loopLink=node->next.first==next? node->next.second : node->prev.second;
790 const auto segment=loopLink->networkLink();
791 if(segment)
792 {
793 const VectorDim orientedB=(segment->source==node->networkNode? 1.0 : -1.0)*segment->burgers();
794 if(segment->isBoundarySegment() || segment->isGrainBoundarySegment() ||
795 (orientedB-link->burgers()).norm()>FLT_EPSILON ||
796 !physicalLinks.insert(segment.get()).second)
797 {
798 valid=false;
799 }
800 if(segment==link)
801 {
802 midpointAfter=k;
803 }
804 }
805 }
806 if(!valid)
807 {
808 break;
809 }
810 }
811 if(!valid || midpointAfter==path.size())
812 {
813 continue;
814 }
815
816 // Commit only after all contributors, translations and target-plane
817 // metadata have been checked. The return side shares the original
818 // network nodes, cancelling the junction on every periodic piece.
819 const auto seed=DN.loops().create(slip->s.cartesian(),glidePlane);
820 std::vector<std::shared_ptr<LoopNodeType>> nodes;
821 for(auto iter=prepared.rbegin();iter!=prepared.rend();++iter)
822 {
823 nodes.push_back(DN.loopNodes().create(seed,iter->node,iter->position,iter->patch,iter->edges));
824 }
825 for(size_t k=0;k+1<prepared.size();++k)
826 {
827 const auto& point=prepared[k];
828 if(k>0)
829 {
830 const auto freeNode=DN.networkNodes().clone(point.node->sID);
831 nodes.push_back(DN.loopNodes().create(seed,freeNode,point.position,point.patch,point.edges));
832 }
833 if(k==midpointAfter)
834 {
835 const auto freeNode=DN.networkNodes().create(midpoint);
836 nodes.push_back(DN.loopNodes().create(seed,freeNode,midpoint,periodicPlane->getPatch(VectorDim::Zero()),EdgePair()));
837 }
838 }
839 DN.insertLoop(seed,nodes);
840 ++formedJunctions;
841 }
842 std::cout << "(" << formedJunctions << " junctions)" << magentaColor << " [" << (std::chrono::duration<double>(std::chrono::system_clock::now() - t0)).count() << " sec]" << defaultColor << std::endl;
843 }
844
845 template <int dim>
846 DislocationJunctionFormation<dim>::DislocationJunctionFormation(DislocationNetwork<dim>& DN_in) :
847 /* init */ DN(DN_in)
848 /* init */,maxJunctionIterations(TextFileParser(DN.ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("maxJunctionIterations",true))
849 /* init */,verboseJunctions(TextFileParser(DN.ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("verboseJunctions",true))
850 /* init */,infiniteLineLength(10000.0)
851 {
852
853 }
854
855 template <int dim>
856 void DislocationJunctionFormation<dim>::initFromFile(const std::string& fileName)
857 {
858 collisionTol=TextFileParser(fileName).readScalar<double>("collisionTol",true);
859 }
860
861 template <int dim>
862 void DislocationJunctionFormation<dim>::formJunctions(const double& dx)
863 {
864 #ifdef _OPENMP
865 // const size_t nThreads = omp_get_max_threads();
866 const size_t nThreads = 1;
867 #else
868 const size_t nThreads = 1;
869 #endif
870
871
872
873 size_t nContracted=1;
874 size_t iterations=0;
875 while(nContracted && iterations<maxJunctionIterations)
876 {
877 std::deque<IntersectionTypeContainerType> intersectionContainer;
878 intersectionContainer.resize(nThreads);
879 findIntersections(intersectionContainer,nThreads);
880
881 nContracted=contractJunctions(intersectionContainer);
882 glissileJunctions(dx);
883 iterations++;
884 }
885
886// breakJunctions();
887
888 }
889
890 template <int dim>
891 void DislocationJunctionFormation<dim>::breakJunctions()
892 {
893 for (const auto &nodeIter : DN.networkNodes())
894 {
895 const auto netNode(nodeIter.second.lock()); // shared_ptr to DislocationNode
896
897 // if (link->isSessile() && link->loopLinks().size() > 1) // a junction
898 if ( netNode->loopNodes().size()==2 // a junction
899 && netNode->neighbors().size()==4
900 )
901 {
902 std::cout<<"Detected 4-point junction at node "<<netNode->sID<<std::endl;
903// (*netNode->loopNodes().rbegin())->reset(netNode->clone());
904 }
905 }
906 }
907
908
909 template <int dim>
910 double DislocationJunctionFormation<dim>::collisionTol=10.0;
911
912 template class DislocationJunctionFormation<3>;
913
914}
915#endif
model::DislocationJunctionFormation::DislocationJunctionFormation
DislocationJunctionFormation(DislocationNetwork< dim > &DN_in)
Definition DislocationJunctionFormation.cpp:846
model::DislocationJunctionFormation::maxJunctionIterations
const size_t maxJunctionIterations
The tolerance (in units of distance) used for collision detection.
Definition DislocationJunctionFormation.h:63
model::DislocationJunctionFormation::initFromFile
static void initFromFile(const std::string &fileName)
Definition DislocationJunctionFormation.cpp:856
model::TextFileParser
Definition TextFileParser.h:123
model::TextFileParser::readScalar
Scalar readScalar(const std::string &key, const bool &verbose=false)
Definition TextFileParser.h:238
t
final point of the source segment t
Definition infiniteDisplacement.m:6
a
a
Definition int_R.m:14
b
b
Definition int_R.m:16
x
x
Definition int_R.m:12
create_polycrystal_figure.source
str source
Definition create_polycrystal_figure.py:11
create_polycrystal_figure.s
str s
Definition create_polycrystal_figure.py:12
create_polycrystal_figure.shift
int shift
Definition create_polycrystal_figure.py:32
model::count
size_t InstanceCount< Derived >::count
Definition InstanceCount.h:70
model::collisionTol
double DislocationJunctionFormation< dim >::collisionTol
Definition DislocationJunctionFormation.cpp:910
model::magentaColor
static std::string magentaColor
Definition TerminalColors.h:31
model::defaultColor
static std::string defaultColor
Definition TerminalColors.h:17
model::GlidePlaneKey
LatticePlaneKey< dim > GlidePlaneKey
Definition GlidePlaneKey.h:20
VectorDim
Eigen::Matrix< double, 3, 1 > VectorDim
Definition pyMoDELib.cxx:66
model::SweepPlane::computeIntersectionPairs
void computeIntersectionPairs()
Definition SweepPlane.h:72
edges
m edges
Definition triangle.cpp:15849
dx
REAL dx
Definition triangle.cpp:6539
t0
REAL t0
Definition triangle.cpp:5171
u
REAL u[4]
Definition triangle.cpp:5168
left
left
Definition triangle.cpp:9232
t1
INEXACT REAL t1
Definition triangle.cpp:5170
k
int k
Definition triangle.cpp:3300
success
enum insertvertexresult success
Definition triangle.cpp:8236
first
vertex first
Definition triangle.cpp:8231
current
int current
Definition triangle.cpp:4591
Generated by doxygen 1.15.0