Mechanics of Defect Evolution Library
Loading...
Searching...
No Matches
DislocationNetwork.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_DislocationNetwork_CPP_
9#define model_DislocationNetwork_CPP_
10
11#include <DislocationDynamics/DislocationDynamicsBase.h>
12#include <PeriodicContourTopology.h>
13#include <numbers>
14
15#include <DislocationNetwork.h>
16#include <DislocationNucleation.h>
17#include <algorithm>
18
19
20namespace model
21{
22
23
24 template <int dim>
25 DislocationNetwork<dim>::DislocationNetwork(MicrostructureContainerType& mc) :
26 /* init */ MicrostructureBase<dim>("DislocationDynamics",mc)
27 /* init */,glideStepsSinceLastClimb(0)
28 /* init */,ddBase(this->microstructures.ddBase)
29 /* init */,networkRemesher(*this)
30 /* init */,junctionsMaker(*this)
31 /* init */,crossSlipModel(DislocationCrossSlip<dim>::getModel(ddBase.poly,ddBase.simulationParameters.traitsIO))
32 /* init */,crossSlipMaker(*this)
33 /* init */,nodeContractor(*this)
34 /* init */,timeStepper(*this)
35 /* init */,stochasticForceGenerator(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("useStochasticForce",true)? new StochasticForceGenerator(ddBase.simulationParameters.traitsIO) : nullptr)
36 /* init */,bulkNucleationModel(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("bulkNucleationModel",true))
37 /* init */,surfaceNucleationModel(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("surfaceNucleationModel",true))
38 /* init */,computeDDinteractions(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("computeDDinteractions",true))
39 /* init */,outputQuadraturePoints(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputQuadraturePoints",true))
40 /* init */,outputLinkingNumbers(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputLinkingNumbers",true))
41 /* init */,outputLoopLength(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputLoopLength",true))
42 /* init */,outputSegmentPairDistances(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputSegmentPairDistances",true))
43 /* init */,outputPlasticDistortionPerSlipSystem(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputPlasticDistortionPerSlipSystem",true))
44 /* init */,outputDislocationDensityPerSlipSystem(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("outputDislocationDensityPerSlipSystem",true))
45 /* init */,computeElasticEnergyPerLength(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("computeElasticEnergyPerLength",true))
46 /* init */,alphaLineTension(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<double>("alphaLineTension",true))
47 /* init */,use_velocityFilter(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<double>("use_velocityFilter",true))
48 /* init */,velocityReductionFactor(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<double>("velocityReductionFactor",true))
49 /* init */,nodalVelocityConstraints(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readMatrixCols<double>("nodalVelocityConstraints",dim,true))
50 /* init */,subcyclingBins(getSubCyclingSet(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readArray<int>("subcyclingBins",true)))
51 /* init */,verboseDislocationNode(TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("verboseDislocationNode",true))
52 {
53 assert(velocityReductionFactor>0.0 && velocityReductionFactor<=1.0);
54 LoopNetworkType::verboseLevel=TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("verboseLoopNetwork",true);
55 verboseDislocationNetwork=TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readScalar<int>("verboseDislocationNetwork",true);
56 }
57
58template <int dim>
59std::set<int> DislocationNetwork<dim>::getSubCyclingSet(const std::vector<int> &inpVector)
60{
61 std::set<int> temp;
62 for (const auto &iv : inpVector)
63 {
64 temp.insert(iv);
65 }
66 return (temp.size()>1)? temp : std::set<int>();
67}
68
69 template <int dim>
70 void DislocationNetwork<dim>::initializeConfiguration(const DefectiveCrystalConfigIO<dim>& configIO,const std::ofstream&,const std::ofstream&)
71 {
72 this->lastUpdateTime=this->microstructures.ddBase.simulationParameters.totalTime;
73
74 LoopType::initFromFile(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt");
75 LoopNodeType::initFromFile(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt");
76 LoopLinkType::initFromFile(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt");
77 NetworkLinkType::initFromFile(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt");
78 DislocationFieldBase<dim>::initFromFile(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt");
79
80 glideSolver=DislocationGlideSolverFactory<dim>::getGlideSolver(*this,TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readString("glideSolverType",false));
81 climbSolver=DislocationClimbSolverFactory<dim>::getClimbSolver(*this,TextFileParser(ddBase.simulationParameters.traitsIO.inputFilesFolder+"/DD.txt").readString("climbSolverType",false));
82 _inclusions=this->microstructures.template getUniqueTypedMicrostructure<InclusionMicrostructure<dim>>();
83
84 setConfiguration(configIO);
85 }
86
87 template <int dim>
88 void DislocationNetwork<dim>::addConfiguration(const DefectiveCrystalConfigIO<dim>& evl)
89 {
90 // DislocationNode<3>::force_count(0);
91 // DislocationLoopNode<3>::force_count(0);
92 // DislocationLoop<3>::force_count(0);
93 // EshelbyInclusionBase<3>::force_count(0);
94 //
95 // this->loopLinks().clear(); // erase base network to clear current config
96
97 const size_t initialLoopCount(DislocationLoop<3>::get_count());
98 const size_t initialNetworkNodeCount(NetworkNodeType::get_count());
99 const size_t initialLoopNodeCount(DislocationLoopNode<3>::get_count());
100
101 // Create Loops
102 std::deque<std::shared_ptr<LoopType>> tempLoops; // keep loops alive during setConfiguration
103 size_t loopNumber=1;
104 for(const auto& loop : evl.loops())
105 {
106 const bool faulted(ddBase.poly.grain(loop.grainID)->rationalLatticeDirection(loop.B).rat.asDouble()!=1.0? true : false);
107 VerboseDislocationNetwork(1,"Creating DislocationLoop "<<loop.sID<<" ("<<loopNumber<<" of "<<evl.loops().size()<<"), type="<<loop.loopType<<", faulted="<<faulted<<", |b|="<<loop.B.norm()<<std::endl;);
108
109 //std::cout<<"Creating DislocationLoop "<<loop.sID<<" ("<<loopNumber<<" of "<<evl.loops().size()<<"), type="<<loop.loopType<<", faulted="<<faulted<<", |b|="<<loop.B.norm()<<std::endl;
110 const size_t loopIDtoUse(loop.sID+initialLoopCount);
111 LoopType::set_count(loopIDtoUse);
112
113 GlidePlaneKey<dim> loopPlaneKey(loop.P, ddBase.poly.grain(loop.grainID)->reciprocalLatticeDirection(loop.N));
114 tempLoops.push_back(this->loops().create(loop.B, ddBase.dislocationDynamics().glidePlaneFactory.getFromKey(loopPlaneKey)));
115 assert(this->loops().get(loopIDtoUse)->sID == loopIDtoUse);
116 loopNumber++;
117 }
118
119 // Create NetworkNodes
120 std::deque<std::shared_ptr<NetworkNodeType>> tempNetNodes; // keep loops alive during setConfiguration
121 size_t netNodeNumber=1;
122 for(const auto& node : evl.nodes())
123 {
124 VerboseDislocationNetwork(1,"Creating DislocationNode "<<node.sID<<" ("<<netNodeNumber<<" of "<<evl.nodes().size()<<")"<<std::endl;);
125 const size_t nodeIDtoUse(node.sID+initialNetworkNodeCount);
126 NetworkNodeType::set_count(nodeIDtoUse);
127 tempNetNodes.push_back(this->networkNodes().create(node.P,node.V,node.climbVelocityScalar,node.velocityReduction));
128 assert(this->networkNodes().get(nodeIDtoUse)->sID==nodeIDtoUse);
129 netNodeNumber++;
130 }
131
132 // Create LoopNodes
133 std::deque<std::shared_ptr<LoopNodeType>> tempLoopNodes; // keep loops alive during setConfiguration
134 size_t loopNodeNumber=1;
135 for(const auto& node : evl.loopNodes())
136 {
137 VerboseDislocationNetwork(1,"Creating DislocationLoopNode "<<node.sID<<" ("<<loopNodeNumber<<" of "<<evl.loopNodes().size()<<")"<<std::endl;);
138 const size_t nodeIDtoUse(node.sID+initialLoopNodeCount);
139 LoopNodeType::set_count(nodeIDtoUse);
140 const auto loop(this->loops().get(node.loopID+initialLoopCount));
141 const auto netNode(this->networkNodes().get(node.networkNodeID+initialNetworkNodeCount));
142 const auto periodicPatch(loop->periodicGlidePlane? loop->periodicGlidePlane->patches().getFromKey(node.periodicShift) : nullptr);
143 std::pair<std::shared_ptr<PeriodicPlaneEdge<dim>>,std::shared_ptr<PeriodicPlaneEdge<dim>>> periodicPatchEdge;
144 if(periodicPatch && node.edgeIDs.first>=0)
145 {
146 const VectorDim position=node.P+node.periodicShift;
147 const auto& edges=periodicPatch->edges();
148 // Old snapshots can have a small loop/network position mismatch.
149 // Use that recorded discrepancy as the reconstruction uncertainty.
150 const double tolerance=FLT_EPSILON+(position-netNode->get_P()).norm();
151 const auto matches=[&](int index)
152 {
153 return index>=0 && size_t(index)<edges.size() && edges[index]->meshIntersection->distanceTo(position)<=tolerance;
154 };
155 if(matches(node.edgeIDs.first) && (node.edgeIDs.second<0 || (node.edgeIDs.second!=node.edgeIDs.first && matches(node.edgeIDs.second))))
156 {
157 periodicPatchEdge.first=edges[node.edgeIDs.first];
158 if(node.edgeIDs.second>=0)
159 {
160 periodicPatchEdge.second=edges[node.edgeIDs.second];
161 }
162 }
163 else
164 {
165 // Older configurations stored contour-vector indices. Those
166 // indices change when its starting edge changes. Recover the
167 // incidence from the saved position, never snap to an unrelated
168 // indexed edge. This is IO compatibility, not contour ordering.
169 std::vector<std::shared_ptr<PeriodicPlaneEdge<dim>>> candidates;
170 for(const auto& edge : edges)
171 {
172 if(edge->meshIntersection->distanceTo(position)<=tolerance)
173 {
174 candidates.push_back(edge);
175 }
176 }
177 const size_t expected=node.edgeIDs.second>=0 ? 2 : 1;
178 if(candidates.size()!=expected)
179 {
180 throw std::runtime_error("Cannot recover saved boundary incidences for loop node "+std::to_string(node.sID));
181 }
182 periodicPatchEdge.first=candidates.front();
183 if(expected==2)
184 {
185 periodicPatchEdge.second=candidates.back();
186 }
187 }
188 }
189 tempLoopNodes.push_back(this->loopNodes().create(loop,netNode,node.P,periodicPatch,periodicPatchEdge));
190 assert(this->loopNodes().get(nodeIDtoUse)->sID==nodeIDtoUse);
191 loopNodeNumber++;
192 }
193
194 // Insert Loops
195 std::map<size_t,std::map<size_t,size_t>> loopMap;
196 for(const auto& looplink : evl.loopLinks())
197 {// Collect LoopLinks by loop IDs
198 loopMap[looplink.loopID].emplace(looplink.sourceID,looplink.sinkID);
199 }
200 assert(loopMap.size()==evl.loops().size());
201
202 for(const auto& loop : evl.loops())
203 {// for each loop in the DefectiveCrystalConfigIO<dim> object
204
205 const auto loopFound=loopMap.find(loop.sID); // there must be an entry with key loopID in loopMap
206 assert(loopFound!=loopMap.end());
207 std::vector<std::shared_ptr<LoopNodeType>> loopNodes;
208 loopNodes.push_back(this->loopNodes().get(loopFound->second.begin()->first+initialLoopNodeCount));
209 for(size_t k=0;k<loopFound->second.size();++k)
210 {
211 const auto nodeFound=loopFound->second.find(loopNodes.back()->sID);
212 if(k<loopFound->second.size()-1)
213 {
214 loopNodes.push_back(this->loopNodes().get(nodeFound->second+initialLoopNodeCount));
215 }
216 else
217 {
218 assert(nodeFound->second+initialLoopNodeCount==loopNodes[0]->sID);
219 }
220 }
221 // std::cout<<" Inserting loop "<<loop.sID<<std::endl;
222 this->insertLoop(this->loops().get(loop.sID+initialLoopCount),loopNodes);
223 // Recover contour bends from older files without changing the IO
224 // format. Geometry is still committed at this point.
225 retainPeriodicContourBends(loopNodes,1.e-8*(ddBase.mesh.xMax()-ddBase.mesh.xMin()).norm());
226 }
227 updateGeometry();
228 }
229
230 template <int dim>
231 void DislocationNetwork<dim>::setConfiguration(const DefectiveCrystalConfigIO<dim>& evl)
232 {
233 DislocationNode<3>::force_count(0);
234 DislocationLoopNode<3>::force_count(0);
235 DislocationLoop<3>::force_count(0);
236 // EshelbyInclusionBase<3>::force_count(0);
237
238 this->loopLinks().clear(); // erase base network to clear current config
239 addConfiguration(evl);
240
241 // // Create Loops
242 // std::deque<std::shared_ptr<LoopType>> tempLoops; // keep loops alive during setConfiguration
243 // size_t loopNumber=1;
244 // for(const auto& loop : evl.loops())
245 // {
246 // const bool faulted(ddBase.poly.grain(loop.grainID)->rationalLatticeDirection(loop.B).rat.asDouble()!=1.0? true : false);
247 // VerboseDislocationNetwork(1,"Creating DislocationLoop "<<loop.sID<<" ("<<loopNumber<<" of "<<evl.loops().size()<<"), type="<<loop.loopType<<", faulted="<<faulted<<", |b|="<<loop.B.norm()<<std::endl;);
248 //
249 // //std::cout<<"Creating DislocationLoop "<<loop.sID<<" ("<<loopNumber<<" of "<<evl.loops().size()<<"), type="<<loop.loopType<<", faulted="<<faulted<<", |b|="<<loop.B.norm()<<std::endl;
250 // const size_t loopIDinFile(loop.sID);
251 // LoopType::set_count(loopIDinFile);
252 //
253 // GlidePlaneKey<dim> loopPlaneKey(loop.P, ddBase.poly.grain(loop.grainID)->reciprocalLatticeDirection(loop.N));
254 // tempLoops.push_back(this->loops().create(loop.B, ddBase.dislocationDynamics().glidePlaneFactory.getFromKey(loopPlaneKey)));
255 // assert(this->loops().get(loopIDinFile)->sID == loopIDinFile);
256 // loopNumber++;
257 // }
258 //
259 // // Create NetworkNodes
260 // std::deque<std::shared_ptr<NetworkNodeType>> tempNetNodes; // keep loops alive during setConfiguration
261 // size_t netNodeNumber=1;
262 // for(const auto& node : evl.nodes())
263 // {
264 // VerboseDislocationNetwork(1,"Creating DislocationNode "<<node.sID<<" ("<<netNodeNumber<<" of "<<evl.nodes().size()<<")"<<std::endl;);
265 // const size_t nodeIDinFile(node.sID);
266 // NetworkNodeType::set_count(nodeIDinFile);
267 // tempNetNodes.push_back(this->networkNodes().create(node.P,node.V,node.climbVelocityScalar,node.velocityReduction));
268 // assert(this->networkNodes().get(nodeIDinFile)->sID==nodeIDinFile);
269 // netNodeNumber++;
270 // }
271 //
272 // // Create LoopNodes
273 // std::deque<std::shared_ptr<LoopNodeType>> tempLoopNodes; // keep loops alive during setConfiguration
274 // size_t loopNodeNumber=1;
275 // for(const auto& node : evl.loopNodes())
276 // {
277 // VerboseDislocationNetwork(1,"Creating DislocationLoopNode "<<node.sID<<" ("<<loopNodeNumber<<" of "<<evl.loopNodes().size()<<")"<<std::endl;);
278 // const size_t nodeIDinFile(node.sID);
279 // LoopNodeType::set_count(nodeIDinFile);
280 // const auto loop(this->loops().get(node.loopID));
281 // const auto netNode(this->networkNodes().get(node.networkNodeID));
282 // const auto periodicPatch(loop->periodicGlidePlane? loop->periodicGlidePlane->patches().getFromKey(node.periodicShift) : nullptr);
283 // const auto periodicPatchEdge((periodicPatch && node.edgeIDs.first>=0)? (node.edgeIDs.second>=0 ? std::make_pair(periodicPatch->edges()[node.edgeIDs.first],
284 // periodicPatch->edges()[node.edgeIDs.second]):
285 // std::make_pair(periodicPatch->edges()[node.edgeIDs.first],nullptr)):std::make_pair(nullptr,nullptr));
286 // tempLoopNodes.push_back(this->loopNodes().create(loop,netNode,node.P,periodicPatch,periodicPatchEdge));
287 // assert(this->loopNodes().get(nodeIDinFile)->sID==nodeIDinFile);
288 // loopNodeNumber++;
289 // }
290 //
291 // // Insert Loops
292 // std::map<size_t,std::map<size_t,size_t>> loopMap;
293 // for(const auto& looplink : evl.loopLinks())
294 // {// Collect LoopLinks by loop IDs
295 // loopMap[looplink.loopID].emplace(looplink.sourceID,looplink.sinkID);
296 // }
297 // assert(loopMap.size()==evl.loops().size());
298 //
299 // for(const auto& loop : evl.loops())
300 // {// for each loop in the DefectiveCrystalConfigIO<dim> object
301 //
302 // const auto loopFound=loopMap.find(loop.sID); // there must be an entry with key loopID in loopMap
303 // assert(loopFound!=loopMap.end());
304 // std::vector<std::shared_ptr<LoopNodeType>> loopNodes;
305 // loopNodes.push_back(this->loopNodes().get(loopFound->second.begin()->first));
306 // for(size_t k=0;k<loopFound->second.size();++k)
307 // {
308 // const auto nodeFound=loopFound->second.find(loopNodes.back()->sID);
309 // if(k<loopFound->second.size()-1)
310 // {
311 // loopNodes.push_back(this->loopNodes().get(nodeFound->second));
312 // }
313 // else
314 // {
315 // assert(nodeFound->second==loopNodes[0]->sID);
316 // }
317 // }
318 // // std::cout<<" Inserting loop "<<loop.sID<<std::endl;
319 // this->insertLoop(this->loops().get(loop.sID),loopNodes);
320 // }
321 // updateGeometry();
322 }
323
324 template <int dim>
325 void DislocationNetwork<dim>::updateGeometry()
326 {
327 VerboseDislocationNetwork(2,"DislocationNetwork::updateGeometry"<<std::endl;);
328 for(auto& loop : this->loops())
329 {// copmute slipped areas and right-handed normal // TODO: PARALLELIZE THIS LOOP
330 loop.second.lock()->updateGeometry();
331 }
332 // updatePlasticDistortionRateFromAreas();
333 VerboseDislocationNetwork(3,"DislocationNetwork::updateGeometry DONE"<<std::endl;);
334 }
335
336 template <int dim>
337 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::averagePlasticDistortion() const
338 {
339 MatrixDim temp(MatrixDim::Zero());
340 for(const auto& loop : this->loops())
341 {
342 temp+= loop.second.lock()->averagePlasticDistortion();
343 }
344 return temp;
345 }
346
347 template <int dim>
348 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::averagePlasticStrain() const
349 {
351 const MatrixDim apd(averagePlasticDistortion());
352 return 0.5*(apd+apd.transpose());
353 }
354
355 template <int dim>
356 std::map<std::pair<int,int>,double> DislocationNetwork<dim>::slipSystemAveragePlasticDistortion() const
357 {
358 std::map<std::pair<int,int>,double> temp; // <grainID,slipSystemID>
359 for(const auto& grain : ddBase.poly.grains)
360 {
361 for(const auto& slipSystem : grain.second->slipSystems())
362 {
363 const std::pair<int,int> key(std::make_pair(grain.first,slipSystem.second->sID));
364 temp[key]=0.0;
365 }
366 }
367
368 for(const auto& weakloop : this->loops())
369 {
370 const auto loop(weakloop.second.lock());
371 if(loop->slipSystem())
372 {
373 const double val = loop->slipSystem()->unitSlip.transpose()*loop->averagePlasticDistortion()*loop->slipSystem()->unitNormal;
374 const std::pair<int,int> key(std::make_pair(loop->grain.grainID,loop->slipSystem()->sID));
375 auto iter(temp.find(key));
376 if(iter!=temp.end())
377 {// iter found
378 iter->second += val ;
379 }
380 else
381 {// iter not found
382 temp.emplace(key,val);
383 }
384 }
385 }
386 return temp;
387 }
388
389 template <int dim>
390 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::averagePlasticDistortionRate() const
391 {
392 MatrixDim temp(MatrixDim::Zero());
393 for(const auto& loop : this->loops())
394 {
395 temp+= loop.second.lock()->averagePlasticDistortionRate();
396 }
397 return temp;
398 }
399
400 template <int dim>
401 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::averagePlasticStrainRate() const
402 {
404 const MatrixDim apdr(averagePlasticDistortionRate());
405 return 0.5*(apdr+apdr.transpose());
406 }
407
408 template <int dim>
409 std::vector<std::tuple<double,double,double,double>> DislocationNetwork<dim>::networkLengthPerSlipSystem() const
410 {
411 std::vector<std::tuple<double,double,double,double>> temp;
412
413 if(ddBase.poly.grains.size())
414 {
415 temp.resize(ddBase.poly.grains.begin()->second->slipSystems().size(),std::make_tuple(0.0,0.0,0.0,0.0));
416 // std::vector<std::tuple<double,double,double,double>> temp(ddBase.poly.grains.begin()->second->slipSystems().size(),std::make_tuple(0.0,0.0,0.0,0.0));
417 for(auto& loop : this->loops())
418 {
419 if(loop.second.lock()->slipSystem())
420 {
421 const size_t ssID(loop.second.lock()->slipSystem()->sID);
422 double& bulkGlissileLength(std::get<0>(temp[ssID]));
423 double& bulkSessileLength(std::get<1>(temp[ssID]));
424 double& boundaryLength(std::get<2>(temp[ssID]));
425 double& grainBoundaryLength(std::get<3>(temp[ssID]));
426
427 for(const auto& loopLink : loop.second.lock()->loopLinks())
428 {
429 if(loopLink->networkLink())
430 {
431 if(!loopLink->networkLink()->hasZeroBurgers())
432 {
433 if(loopLink->networkLink()->isGrainBoundarySegment())
434 {
435 grainBoundaryLength+=loopLink->networkLink()->chord().norm();
436 }
437 else if(loopLink->networkLink()->isBoundarySegment())
438 {
439 boundaryLength+=loopLink->networkLink()->chord().norm();
440 }
441 else
442 {
443 if(loopLink->networkLink()->isSessile())
444 {
445 bulkSessileLength+=loopLink->networkLink()->chord().norm()/loopLink->networkLink()->loopLinks().size();
446 }
447 else
448 {
449 bulkGlissileLength+=loopLink->networkLink()->chord().norm()/loopLink->networkLink()->loopLinks().size();
450 }
451 }
452 }
453 }
454 }
455 }
456 }
457 }
458 return temp;
459 }
460
461 template <int dim>
462 std::tuple<double,double,double,double> DislocationNetwork<dim>::networkLength() const
463 {
469 double bulkGlissileLength(0.0);
470 double bulkSessileLength(0.0);
471 double boundaryLength(0.0);
472 double grainBoundaryLength(0.0);
473
474 for(auto& loop : this->loops())
475 {
476 for(const auto& loopLink : loop.second.lock()->loopLinks())
477 {
478 if(loopLink->networkLink())
479 {
480 if(!loopLink->networkLink()->hasZeroBurgers())
481 {
482 if(loopLink->networkLink()->isGrainBoundarySegment())
483 {
484 grainBoundaryLength+=loopLink->networkLink()->chord().norm();
485 }
486 else if(loopLink->networkLink()->isBoundarySegment())
487 {
488 boundaryLength+=loopLink->networkLink()->chord().norm();
489 }
490 else
491 {
492 if(loopLink->networkLink()->isSessile())
493 {
494 bulkSessileLength+=loopLink->networkLink()->chord().norm()/loopLink->networkLink()->loopLinks().size();
495 }
496 else
497 {
498 bulkGlissileLength+=loopLink->networkLink()->chord().norm()/loopLink->networkLink()->loopLinks().size();
499 }
500 }
501 }
502 }
503 }
504 }
505 return std::make_tuple(bulkGlissileLength,bulkSessileLength,boundaryLength,grainBoundaryLength);
506 }
507
508
509 template <int dim>
510 const std::shared_ptr<InclusionMicrostructure<dim>>& DislocationNetwork<dim>::inclusions() const
511 {
512 return _inclusions;
513 }
514
515 template <int dim>
516 bool DislocationNetwork<dim>::contract(std::shared_ptr<NetworkNodeType> nA,
517 std::shared_ptr<NetworkNodeType> nB)
518 {
519 return nodeContractor.contract(nA,nB);
520 }
521
522 template <int dim>
523 typename DislocationNetwork<dim>::VectorDim DislocationNetwork<dim>::displacement(const VectorDim& x,const NodeType* const,const ElementType* const,const SimplexDim* const) const
524 {
529 VectorDim temp(VectorDim::Zero());
530
531 for(const auto& loop : this->loops())
532 {// sum solid angle of each loop
533 for(const auto& shift : ddBase.periodicShifts)
534 {
535 temp-=loop.second.lock()->solidAngle(x+shift)/4.0/std::numbers::pi*loop.second.lock()->burgers();
536 }
537 }
538
539 for(const auto& link : this->networkLinks())
540 {// sum line-integral part of displacement field per segment
541 if( !link.second.lock()->hasZeroBurgers())
542 {
543 for(const auto& shift : ddBase.periodicShifts)
544 {
545 temp+=link.second.lock()->straight.displacement(x+shift);
546 }
547 }
548 }
549
550 return temp;
551 }
552
553 template <int dim>
554 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::averageStress() const
555 {
560 return MatrixDim::Zero();
561 }
562
563 template <int dim>
564 typename DislocationNetwork<dim>::VectorMSize DislocationNetwork<dim>::mobileConcentration(const VectorDim& x, const NodeType* const node, const ElementType* const ele,const SimplexDim* const guess) const
565 {
566 VectorMSize temp(VectorMSize::Zero());
567 if(climbSolver)
568 {
569 const auto pointGrains(this->pointGrains(x,node,ele,guess));
570 // std::set<const Grain<dim>*> pointGrains;
571 // if(node)
572 // {
573 // for(const auto& nodeEle : *node)
574 // {
575 // pointGrains.emplace(&ddBase.poly.grain(nodeEle->simplex.region->regionID));
576 // }
577 // }
578 // else
579 // {
580 // if(ele)
581 // {
582 // pointGrains.emplace(&ddBase.poly.grain(ele->simplex.region->regionID));
583 // }
584 // else
585 // {
586 // const std::pair<bool,const Simplex<dim,dim>*> found(ddBase.mesh.searchWithGuess(x,guess));
587 // if(found.first)
588 // {
589 // pointGrains.emplace(&ddBase.poly.grain(found.second->region->regionID));
590 // }
591 // }
592 // }
593
594 if(pointGrains.size())
595 {
596 for(const auto& link : this->networkLinks())
597 {// sum stress field per segment
598 if( !link.second.lock()->hasZeroBurgers()
599 )
600 {
601
602 const auto segment(link.second.lock());
603 const auto segGrains(segment->grains());
604 std::set<const Grain<dim>*> intersect;
605 std::set_intersection(pointGrains.begin(), pointGrains.end(), segGrains.begin(), segGrains.end(),
606 std::inserter(intersect, intersect.begin()));
607 if(intersect.size()==1)
608 {
609 // const int grainID((*intersect.begin())->region.regionID);
610 // StressStraight<3> ss(ddBase.poly,segment->source->get_P(),segment->sink->get_P(),segment->burgers(),ddBase.dislocationDynamics().EwaldLength);
611 // for(const auto& shift : ddBase.periodicShifts)
612 // {
613 temp+=segment->clusterConcentration(x,climbSolver->CD->cdp);
614
615 // temp+=ss.clusterConcentration(x+shift,grainID, segment->source->climbDirection(), segment->source->climbVelocityScalar , segment->sink->climbDirection(), segment->sink->climbVelocityScalar, climbSolver->CD->cdp);
616 // }
617 }
618 }
619 }
620 }
621 }
622 return temp;
623 }
624
625
626 template<int dim>
627 typename DislocationNetwork<dim>::VectorISize DislocationNetwork<dim>::immobileClusters(const VectorDim&,const NodeType* const,const ElementType* const,const SimplexDim* const) const
628 {
629 return VectorISize::Zero();
630 }
631
632
633 template <int dim>
634 typename DislocationNetwork<dim>::MatrixDim DislocationNetwork<dim>::stress(const VectorDim& x,const NodeType* const,const ElementType* const,const SimplexDim* const) const
635 {
640 MatrixDim temp(MatrixDim::Zero());
641 for(const auto& link : this->networkLinks())
642 {// sum stress field per segment
643 if( !link.second.lock()->hasZeroBurgers()
644 )
645 {
646 for(const auto& shift : ddBase.periodicShifts)
647 {
648 temp+=link.second.lock()->straight.stress(x+shift);
649 }
650 }
651 }
652 return temp;
653 }
654
655 template <int dim>
656 typename DislocationNetwork<dim>::VectorDim DislocationNetwork<dim>::inelasticDisplacementRate(const VectorDim&, const NodeType* const, const ElementType* const,const SimplexDim* const) const
657 {
658 return VectorDim::Zero();
659 }
660
661 template <int dim>
662 void DislocationNetwork<dim>::updateConfiguration()
663 {
664 this->lastUpdateTime=this->microstructures.ddBase.simulationParameters.totalTime;
665
666 moveNodes(ddBase.simulationParameters.dt);
667 executeSingleGlideStepDiscreteEvents(ddBase.simulationParameters.runID);
668 if(isClimbStep())
669 {
670 glideStepsSinceLastClimb=0;
671 }
672 else
673 {
674 glideStepsSinceLastClimb++;
675 }
676 updateGeometry();
677 }
678
679 template <int dim>
680 double DislocationNetwork<dim>::getDt() const
681 {
682 if(isClimbStep())
683 {
684 return timeStepper.getDt(climbSolver->vClimbRef,1.0e-7);
685
686 }
687 else
688 {
689 // climbSolver->CD->cdp
690 return timeStepper.getDt(ddBase.poly.cs,1.0e-7);
691 }
692 }
693
694 template <int dim>
695 bool DislocationNetwork<dim>::isClimbStep() const
696 {
697 if(glideSolver)
698 {
699 return climbSolver? (glideStepsSinceLastClimb && averagePlasticDistortionRate().norm()<climbSolver->glideEquilibriumRate): false;
700 }
701 else
702 {
703 if(climbSolver)
704 {
705 return true;
706 }
707 else
708 {
709 return false;
710 }
711 }
712 }
713
714template <int dim>
715void DislocationNetwork<dim>::reSolve()
716{
717
718 const bool isClimbingStep(isClimbStep());
719 if(isClimbingStep)
720 {
721 std::cout<<" climbStep"<<std::flush;
722 }
723 else
724 {
725 std::cout<<" glideStep"<<std::flush;
726 }
727 double maxVelocity = 0.0;
728 for (const auto &nodeIter : this->networkNodes())
729 {
730 const double vNorm(nodeIter.second.lock()->get_V().norm());
731 if (vNorm > maxVelocity)
732 {
733 maxVelocity = vNorm;
734 }
735 }
736
737#ifdef _OPENMP
738 const size_t nThreads = omp_get_max_threads();
739#else
740 const size_t nThreads = 1;
741#endif
742
744 std::map<int, int> velocityBinMap{{1,0}};
745 for (const auto &binVal : subcyclingBins)
746 {
747 velocityBinMap.emplace(binVal, 0);
748 }
749
750
751 std::cout <<" ,updating updateForcesAndVelocities (" << nThreads << " threads) " << std::flush;
752#ifdef _OPENMP
753 {
754 // Topology is fixed during field evaluation. Preserve container order
755 // and shared ownership, but make OpenMP's indexed access constant-time.
756 // Release the snapshot before any subsequent topology changes.
757 std::vector<std::shared_ptr<NetworkLinkType>> fieldLinks;
758 fieldLinks.reserve(this->networkLinks().size());
759 for(const auto& entry : this->networkLinks())
760 {
761 fieldLinks.push_back(entry.second.lock());
762 }
763#pragma omp parallel for
764 for (size_t k = 0; k < fieldLinks.size(); ++k)
765 {
766 const auto& link=fieldLinks[k];
767 const int velGroup((subcyclingBins.size() && !isClimbingStep) ? link->velocityGroup(maxVelocity, subcyclingBins) : 1);
768
769 if ((ddBase.simulationParameters.runID % velGroup) == 0)
770 {
771 for(auto& qp : link->quadraturePoints())
772 {
773 qp.updateForcesAndVelocities(*link,isClimbingStep);
774 }
775 }
776 // else
777 // {
778 // link->assembleGlide(false);
779 // }
780 }
781 }
782
783#else
784 for (auto &linkIter : this->networkLinks())
785 {
786 const auto link=linkIter.second.lock();
787 const int velGroup(subcyclingBins.size() ? linkIter.second.lock()->velocityGroup(maxVelocity, subcyclingBins) : 1);
788
789 if ((ddBase.simulationParameters.runID % velGroup) == 0)
790 {
791 for(auto& qp : link->quadraturePoints())
792 {
793 qp.updateForcesAndVelocities(*link,isClimbingStep);
794 }
795 }
796 // else
797 // {
798 // linkIter.second.lock()->assembleGlide(false);
799 // }
800 }
801#endif
802
803 Eigen::VectorXd X(Eigen::VectorXd::Zero(0));
804 if(isClimbingStep)
805 {
806 climbSolver->computeClimbScalarVelocities(false);
807 size_t k=0;
808 for (auto& networkNode : this->networkNodes())
809 {
810 networkNode.second.lock()->climbVelocityScalar=climbSolver->scalarVelocities()[k];
811 ++k;
812 }
813 X=climbSolver->getNodeVelocities();
814 }
815 else
816 {
817 if(glideSolver)
818 {
819 size_t k=0;
820 for (auto& networkNode : this->networkNodes())
821 {
822 networkNode.second.lock()->climbVelocityScalar.setZero();
823 ++k;
824 }
825 X=glideSolver->getNodeVelocities();
826 }
827 }
828 if(int(NdofXnode*this->networkNodes().size())==X.size())
829 {
830 size_t k=0;
831 for (auto& networkNode : this->networkNodes())
832 {
833 networkNode.second.lock()->set_V(X.segment(NdofXnode*k,NdofXnode),isClimbingStep); // double cast to remove some numerical noise
834 ++k;
835 }
836 }
837 else
838 {
839 std::cout<<"NdofXnode*this->networkNodes().size()="<<NdofXnode*this->networkNodes().size()<<std::endl;
840 std::cout<<"vSolver->getNodeVelocities().size()="<<X.size()<<std::endl;
841 throw std::runtime_error("vSolver returned wrong velocity vector size.");
842 }
843
844 VerboseDislocationNetwork(2,"DislocationNetwork::updateRates"<<std::endl;);
845 for(auto& loop : this->loops())
846 {// copmute slipped areas and right-handed normal // TODO: PARALLELIZE THIS LOOP
847 loop.second.lock()->updateRates();
848 }
849 // updatePlasticDistortionRateFromAreas();
850 VerboseDislocationNetwork(3,"DislocationNetwork::updateRates DONE"<<std::endl;);
851
852// storeSingleGlideStepDiscreteEvents(ddBase.simulationParameters.runID);
853
854}
855
856 template <int dim>
857 void DislocationNetwork<dim>::solve()
858 {
859
860 const bool isClimbingStep(isClimbStep());
861 if(isClimbingStep)
862 {
863 std::cout<<" climbStep"<<std::flush;
864 }
865 else
866 {
867 std::cout<<" glideStep"<<std::flush;
868 }
869 double maxVelocity = 0.0;
870 for (const auto &nodeIter : this->networkNodes())
871 {
872 const double vNorm(nodeIter.second.lock()->get_V().norm());
873 if (vNorm > maxVelocity)
874 {
875 maxVelocity = vNorm;
876 }
877 }
878
879 #ifdef _OPENMP
880 const size_t nThreads = omp_get_max_threads();
881 #else
882 const size_t nThreads = 1;
883 #endif
884
886 std::map<int, int> velocityBinMap{{1,0}};
887 for (const auto &binVal : subcyclingBins)
888 {
889 velocityBinMap.emplace(binVal, 0);
890 }
891
892 std::cout <<" creating qPoints "<< std::flush;
893 for (const auto &links : this->networkLinks())
894 {
895
896 const int velGroup((subcyclingBins.size() && !isClimbingStep) ? links.second.lock()->velocityGroup(maxVelocity, subcyclingBins) : 1);
897 auto velocityBinIter(velocityBinMap.find(velGroup));
898 assert(velocityBinIter != velocityBinMap.end());
899 velocityBinIter->second++;
900
901 if ((ddBase.simulationParameters.runID % velGroup) == 0)
902 {
903 links.second.lock()->createQuadraturePoints(isClimbingStep);
904 }
905 }
906
907 std::cout <<" ,updating qPoints (" << nThreads << " threads) " << std::flush;
908 #ifdef _OPENMP
909 {
910 // Topology is fixed during field evaluation. Preserve container order
911 // and shared ownership, but make OpenMP's indexed access constant-time.
912 // Release the snapshot before any subsequent topology changes.
913 std::vector<std::shared_ptr<NetworkLinkType>> fieldLinks;
914 fieldLinks.reserve(this->networkLinks().size());
915 for(const auto& entry : this->networkLinks())
916 {
917 fieldLinks.push_back(entry.second.lock());
918 }
919 #pragma omp parallel for
920 for (size_t k = 0; k < fieldLinks.size(); ++k)
921 {
922 const auto& link=fieldLinks[k];
923 const int velGroup((subcyclingBins.size() && !isClimbingStep) ? link->velocityGroup(maxVelocity, subcyclingBins) : 1);
924
925 if ((ddBase.simulationParameters.runID % velGroup) == 0)
926 {
927 link->updateQuadraturePoints(isClimbingStep);
928 }
929 // else
930 // {
931 // link->assembleGlide(false);
932 // }
933 }
934 }
935
936 #else
937 for (auto &linkIter : this->networkLinks())
938 {
939 const int velGroup(subcyclingBins.size() ? linkIter.second.lock()->velocityGroup(maxVelocity, subcyclingBins) : 1);
940
941 if ((ddBase.simulationParameters.runID % velGroup) == 0)
942 {
943 linkIter.second.lock()->updateQuadraturePoints(isClimbingStep);
944 }
945 // else
946 // {
947 // linkIter.second.lock()->assembleGlide(false);
948 // }
949 }
950 #endif
951
952 Eigen::VectorXd X(Eigen::VectorXd::Zero(0));
953 if(isClimbingStep)
954 {
955 climbSolver->computeClimbScalarVelocities(true);
956 size_t k=0;
957 for (auto& networkNode : this->networkNodes())
958 {
959 networkNode.second.lock()->climbVelocityScalar=climbSolver->scalarVelocities()[k];
960 ++k;
961 }
962 X=climbSolver->getNodeVelocities();
963 }
964 else
965 {
966 if(glideSolver)
967 {
968 size_t k=0;
969 for (auto& networkNode : this->networkNodes())
970 {
971 networkNode.second.lock()->climbVelocityScalar.setZero();
972 ++k;
973 }
974 X=glideSolver->getNodeVelocities();
975 }
976 }
977 if(int(NdofXnode*this->networkNodes().size())==X.size())
978 {
979 size_t k=0;
980 for (auto& networkNode : this->networkNodes())
981 {
982 networkNode.second.lock()->set_V(X.segment(NdofXnode*k,NdofXnode),isClimbingStep); // double cast to remove some numerical noise
983 ++k;
984 }
985 }
986 else
987 {
988 std::cout<<"NdofXnode*this->networkNodes().size()="<<NdofXnode*this->networkNodes().size()<<std::endl;
989 std::cout<<"vSolver->getNodeVelocities().size()="<<X.size()<<std::endl;
990 throw std::runtime_error("vSolver returned wrong velocity vector size.");
991 }
992
993 VerboseDislocationNetwork(2,"DislocationNetwork::updateRates"<<std::endl;);
994 for(auto& loop : this->loops())
995 {// copmute slipped areas and right-handed normal // TODO: PARALLELIZE THIS LOOP
996 loop.second.lock()->updateRates();
997 }
998 // updatePlasticDistortionRateFromAreas();
999 VerboseDislocationNetwork(3,"DislocationNetwork::updateRates DONE"<<std::endl;);
1000
1001 storeSingleGlideStepDiscreteEvents(ddBase.simulationParameters.runID);
1002
1003 }
1004
1005 template <int dim>
1006 void DislocationNetwork<dim>::moveNodes(const double & dt_in)
1007 {
1009 const auto t0= std::chrono::system_clock::now();
1010 std::cout<<"Moving DislocationNodes (dt="<<dt_in<< ")... "<<std::flush;
1011 // Thermal hazards use this step's accepted time interval, not the
1012 // previous dt still present during the velocity solve.
1013 if(crossSlipModel && crossSlipModel->isThermal() && !isClimbStep())
1014 {
1015 crossSlipMaker.findCrossSlipSegments(dt_in);
1016 crossSlipMaker.execute();
1017 }
1018 std::vector<std::pair<std::shared_ptr<LoopNodeType>,bool>> savedFlags;
1019 for(const auto& entry : this->loopNodes())
1020 {
1021 const auto node=entry.second.lock();
1022 savedFlags.emplace_back(node,node->periodicContourVertex);
1023 }
1024 const auto savedDangling=danglingBoundaryLoopNodes;
1025 for(const auto& entry : this->loops())
1026 {
1027 std::vector<std::shared_ptr<LoopNodeType>> nodes;
1028 for(const auto* link : entry.second.lock()->linkSequence())
1029 {
1030 nodes.push_back(link->source);
1031 }
1032 retainPeriodicContourBends(nodes,1.e-8*(ddBase.mesh.xMax()-ddBase.mesh.xMin()).norm());
1033 }
1034 danglingBoundaryLoopNodes.clear();
1035 std::vector<std::pair<std::shared_ptr<LoopNodeType>,typename LoopNodeType::SplineNodeType>> saved;
1036 for(const auto& entry : this->loopNodes())
1037 {
1038 const auto node=entry.second.lock();
1039 saved.emplace_back(node,static_cast<const typename LoopNodeType::SplineNodeType&>(*node));
1040 }
1041 trialNodeMotion=true;
1042 try
1043 {
1044 // Network positions remain unchanged throughout this loop, so
1045 // shared nodes cannot be advanced twice through different loops.
1046 for(const auto& entry : this->networkNodes())
1047 {
1048 const auto node=entry.second.lock();
1049 node->trySet_P(node->get_P()+node->get_V()*dt_in);
1050 }
1051 updateBoundaryNodes();
1052 trialNodeMotion=false;
1053 }
1054 catch(...)
1055 {
1056 for(const auto& entry : saved)
1057 {
1058 static_cast<typename LoopNodeType::SplineNodeType&>(*entry.first)=entry.second;
1059 }
1060 for(const auto& entry : savedFlags)
1061 {
1062 entry.first->periodicContourVertex=entry.second;
1063 }
1064 danglingBoundaryLoopNodes=savedDangling;
1065 trialNodeMotion=false;
1066 throw;
1067 }
1068 std::cout<<magentaColor<<std::setprecision(3)<<std::scientific<<" ["<<(std::chrono::duration<double>(std::chrono::system_clock::now()-t0)).count()<<" sec]."<<defaultColor<<std::endl;
1069 }
1070
1071 template <int dim>
1072 void DislocationNetwork<dim>::storeSingleGlideStepDiscreteEvents(const long int&)
1073 {
1074 if(!crossSlipModel || !crossSlipModel->isThermal())
1075 {
1076 crossSlipMaker.findCrossSlipSegments();
1077 crossSlipMaker.execute();
1078 }
1079 }
1080
1081 template <int dim>
1082 void DislocationNetwork<dim>::executeSingleGlideStepDiscreteEvents(const long int& runID)
1083 {
1084
1085 // crossSlipMaker.execute();
1086 // //! 13- Node redistribution
1087 networkRemesher.remesh(runID);
1088 // //! 12- Form Junctions
1089 junctionsMaker.formJunctions(3.0*networkRemesher.Lmin);
1090 //Calling remesh again so that any other topological changes created by junctions whihc are otherwise removable can be removed
1091 // //! 13- Node redistribution
1092 // this->io().output(runID);
1093
1094 networkRemesher.remesh(runID);
1095 // updateVirtualBoundaryLoops();
1096
1097 DislocationNucleation<dim> nucleator(*this);
1098 nucleator.bulkNucleate(bulkNucleationModel);
1099 nucleator.surfaceNucleate(surfaceNucleationModel);
1100
1101 }
1102
1103 template <int dim>
1104 void DislocationNetwork<dim>::output(DefectiveCrystalConfigIO<dim>& configIO,DefectiveCrystalAuxIO<dim>& auxIO,std::ofstream& f_file,std::ofstream& F_labels) const
1105 {
1106
1107 for(const auto& loop : this->loops())
1108 {
1109 configIO.loops().emplace_back(*loop.second.lock());
1110 }
1111
1112 // Store LoopNodes
1113 for(const auto& node : this->loopNodes())
1114 {
1115 configIO.loopNodes().emplace_back(*node.second.lock());
1116 }
1117
1118 // Store LoopLinks
1119 for(const auto& link : this->loopLinks())
1120 {
1121 configIO.loopLinks().emplace_back(link.second);
1122 }
1123
1124 // Store NetworkNodes
1125 for(const auto& node : this->networkNodes())
1126 {
1127 configIO.nodes().emplace_back(*node.second.lock());
1128 }
1129
1130 // AuxIO
1131 if (this->outputQuadraturePoints)
1132 {
1133 for (const auto& link : this->networkLinks())
1134 {
1135 for(const auto& qPoint : link.second.lock()->quadraturePoints())
1136 {
1137 auxIO.quadraturePoints().push_back(qPoint);
1138 }
1139 }
1140 }
1141
1142 const std::tuple<double,double,double,double> length(this->networkLength());
1143 const double densityFactor(1.0/this->ddBase.mesh.volume()/std::pow(this->ddBase.poly.b_SI,2));
1144 f_file<<std::get<0>(length)*densityFactor<<" "<<std::get<1>(length)*densityFactor<<" "<<std::get<2>(length)*densityFactor<<" "<<std::get<3>(length)*densityFactor<<" ";
1145 if(ddBase.simulationParameters.runID==0)
1146 {
1147 F_labels<<"glissile density [m^-2]\n";
1148 F_labels<<"sessile density [m^-2]\n";
1149 F_labels<<"boundary density [m^-2]\n";
1150 F_labels<<"grain boundary density [m^-2]\n";
1151 }
1152
1153 if(outputDislocationDensityPerSlipSystem)
1154 {
1155 const auto ssDen(this->networkLengthPerSlipSystem());
1156 for(const auto& tup : ssDen)
1157 {
1158 f_file<<std::get<0>(tup)*densityFactor<<" "<<std::get<1>(tup)*densityFactor<<" "<<std::get<2>(tup)*densityFactor<<" "<<std::get<3>(tup)*densityFactor<<" ";
1159 }
1160
1161 if(ddBase.simulationParameters.runID==0)
1162 {
1163 for(size_t kd=0;kd<ssDen.size();++kd)
1164 {
1165 F_labels<<"SlipSystem_"<<kd<<" glissile density [m^-2]\n";
1166 F_labels<<"SlipSystem_"<<kd<<"sessile density [m^-2]\n";
1167 F_labels<<"SlipSystem_"<<kd<<"boundary density [m^-2]\n";
1168 F_labels<<"SlipSystem_"<<kd<<"grain boundary density [m^-2]\n";
1169 }
1170 }
1171
1172 }
1173
1174 if(this->outputPlasticDistortionPerSlipSystem)
1175 {
1176 const auto ssapd(slipSystemAveragePlasticDistortion());
1177 for(const auto& pair : ssapd)
1178 {
1179 f_file<<pair.second<<" ";
1180 if(ddBase.simulationParameters.runID==0)
1181 {
1182 F_labels<<"betaP_g"+std::to_string(pair.first.first)+"_ss"+std::to_string(pair.first.second)+" [-]\n";
1183 }
1184 }
1185 }
1186
1187 if(this->computeElasticEnergyPerLength)
1188 {
1189 double eE(0.0);
1190 double eC(0.0);
1191
1192 for(const auto& linkIter : this->networkLinks())
1193 {// Collect LoopLinks by loop IDs
1194 const auto link(linkIter.second.lock());
1195 for(const auto& qPoint : link->quadraturePoints())
1196 {
1197 eE+=qPoint.elasticEnergyPerLength*qPoint.dL;
1198 eC+=qPoint.coreEnergyPerLength*qPoint.dL;
1199 }
1200 }
1201 f_file<<eE<<" "<<eC<<" ";
1202 if(ddBase.simulationParameters.runID==0)
1203 {
1204 F_labels<<"dislocation elastic energy [mu b^3]\n";
1205 F_labels<<"dislocation core energy [mu b^3]\n";
1206 }
1207 }
1208 }
1209
1210 template <int dim>
1211 void DislocationNetwork<dim>::updateBoundaryNodes(const std::function<void()>& checkpoint)
1212 {
1213 using NodePtr=std::shared_ptr<LoopNodeType>;
1214 using NetPtr=std::shared_ptr<NetworkNodeType>;
1215 using Edge=PeriodicPlaneEdge<dim>;
1216 using Record=PolygonPatchIntersectionRecord<dim,LoopNodeType>;
1217 const double tolerance=1.e-8*(ddBase.mesh.xMax()-ddBase.mesh.xMin()).norm();
1218 struct Arc
1219 {
1220 NodePtr source,sink;
1221 VectorDim chord;
1222 bool reversed;
1223 std::set<const NetworkLinkType*> links;
1224 std::vector<NodePtr> subdivisions;
1225 size_t group;
1226 };
1227 struct Plan
1228 {
1229 std::shared_ptr<LoopType> loop;
1230 std::vector<NodePtr> oldNodes,anchors,newNodes;
1231 std::vector<size_t> arcs;
1232 std::vector<Record> contour;
1233 };
1234 std::vector<Plan> plans;
1235 std::vector<Arc> arcs;
1236 std::set<NodePtr> promoted;
1237 // Snapshot complete loops first. These strong references also protect
1238 // borrowed provenance pointers until all loops have been reconciled.
1239 for(const auto& entry : this->loops())
1240 {
1241 const auto loop=entry.second.lock();
1242 if(!loop->periodicGlidePlane)
1243 {
1244 continue;
1245 }
1246 Plan plan;
1247 plan.loop=loop;
1248 std::vector<size_t> anchorPositions;
1249 for(const auto* link : loop->linkSequence())
1250 {
1251 plan.oldNodes.push_back(link->source);
1252 // Physical boundary vertices define the retained slipped area;
1253 // only artificial periodic subdivisions may be omitted.
1254 if(!link->source->isPeriodicSubdivision())
1255 {
1256 anchorPositions.push_back(plan.oldNodes.size()-1);
1257 plan.anchors.push_back(link->source);
1258 }
1259 }
1260 if(plan.anchors.size()<3)
1261 {
1262 if(retracesContourSegment(plan.oldNodes,anchorPositions,tolerance))
1263 {
1264 // Disconnect only this loop's memberships. Contraction has
1265 // transported derived points along their original segments;
1266 // both resulting paths have now been verified to retrace.
1267 plans.push_back(std::move(plan));
1268 continue;
1269 }
1270 // No collapse proof: retain the actual contour, including seam
1271 // bends. Never infer zero slipped area from the anchor count.
1272 plan.anchors.clear();
1273 anchorPositions.clear();
1274 for(size_t i=0;i<plan.oldNodes.size();++i)
1275 {
1276 if(plan.anchors.empty() ||
1277 (plan.oldNodes[i]->get_P()-plan.anchors.back()->get_P()).norm()>tolerance)
1278 {
1279 promoted.insert(plan.oldNodes[i]);
1280 plan.anchors.push_back(plan.oldNodes[i]);
1281 anchorPositions.push_back(i);
1282 }
1283 }
1284 if(plan.anchors.size()>1 &&
1285 (plan.anchors.front()->get_P()-plan.anchors.back()->get_P()).norm()<=tolerance)
1286 {
1287 plan.anchors.pop_back();
1288 anchorPositions.pop_back();
1289 }
1290 if(plan.anchors.size()<3)
1291 {
1292 throw std::runtime_error("Cannot establish a retraced periodic contour.");
1293 }
1294 }
1295 for(size_t i=0; i<plan.anchors.size(); ++i)
1296 {
1297 Arc arc;
1298 arc.source=plan.anchors[i];
1299 arc.sink=plan.anchors[(i+1)%plan.anchors.size()];
1300 arc.chord=arc.sink->get_P()-arc.source->get_P();
1301 arc.reversed=arc.source->networkNode->sID>arc.sink->networkNode->sID;
1302 if(arc.source->networkNode==arc.sink->networkNode)
1303 {
1304 for(int axis=0; axis<dim; ++axis)
1305 {
1306 if(std::abs(arc.chord[axis])>tolerance)
1307 {
1308 arc.reversed=arc.chord[axis]<0;
1309 break;
1310 }
1311 }
1312 }
1313 if(arc.reversed)
1314 {
1315 arc.chord=-arc.chord;
1316 }
1317 size_t pos=anchorPositions[i];
1318 const size_t end=anchorPositions[(i+1)%anchorPositions.size()];
1319 do
1320 {
1321 if(const auto link=plan.oldNodes[pos]->next.second->networkLink())
1322 {
1323 arc.links.insert(link.get());
1324 }
1325 pos=(pos+1)%plan.oldNodes.size();
1326 if(pos!=end)
1327 {
1328 arc.subdivisions.push_back(plan.oldNodes[pos]);
1329 }
1330 }
1331 while(pos!=end);
1332 arc.group=arcs.size();
1333 plan.arcs.push_back(arcs.size());
1334 arcs.push_back(std::move(arc));
1335 }
1336 std::vector<std::pair<VectorLowerDim,const LoopNodeType* const>> polygon;
1337 for(const auto& node : plan.anchors)
1338 {
1339 polygon.emplace_back(loop->periodicGlidePlane->referencePlane->localPosition(node->get_P()),node.get());
1340 }
1341 plan.contour=loop->periodicGlidePlane->polygonPatchIntersection(polygon,false);
1342 plans.push_back(std::move(plan));
1343 }
1344 // Establish junction lineage from the pre-edit network, not from
1345 // coincident positions or the endpoints left after trimming.
1346 const auto root=[&](size_t i)
1347 {
1348 while(arcs[i].group!=i)
1349 {
1350 i=arcs[i].group;
1351 }
1352 return i;
1353 };
1354 std::map<std::pair<size_t,size_t>,std::vector<size_t>> endpointArcs;
1355 for(size_t i=0; i<arcs.size(); ++i)
1356 {
1357 const auto key=std::minmax(arcs[i].source->networkNode->sID,arcs[i].sink->networkNode->sID);
1358 auto& peers=endpointArcs[key];
1359 for(size_t j : peers)
1360 {
1361 if((arcs[i].chord-arcs[j].chord).norm()<=tolerance)
1362 {
1363 for(const auto* link : arcs[i].links)
1364 {
1365 if(arcs[j].links.count(link))
1366 {
1367 arcs[root(i)].group=root(j);
1368 break;
1369 }
1370 }
1371 }
1372 }
1373 peers.push_back(i);
1374 }
1375 struct Crossing
1376 {
1377 double parameter;
1378 NetPtr node;
1379 };
1380 std::map<size_t,std::vector<Crossing>> crossings;
1381 for(size_t i=0; i<arcs.size(); ++i)
1382 {
1383 const auto& arc=arcs[i];
1384 const VectorDim d=arc.sink->get_P()-arc.source->get_P();
1385 for(const auto& node : arc.subdivisions)
1386 {
1387 if(!danglingBoundaryLoopNodes.count(node.get()) && d.squaredNorm()>0)
1388 {
1389 double t=d.dot(node->get_P()-arc.source->get_P())/d.squaredNorm();
1390 crossings[root(i)].push_back({arc.reversed ? 1-t : t,node->networkNode});
1391 }
1392 }
1393 }
1394 // Select a representative per original junction before editing any
1395 // positions. Prefer its existing RVE copy when several copies survive.
1396 std::map<size_t,VectorDim> committedPositions;
1397 for(const auto& plan : plans)
1398 {
1399 for(const auto& record : plan.contour)
1400 {
1401 if(const auto index=std::get<6>(record).originalVertex)
1402 {
1403 const auto node=plan.anchors.at(*index)->networkNode;
1404 const VectorDim mapped=plan.loop->periodicGlidePlane->referencePlane->globalPosition(std::get<0>(record))+std::get<1>(record);
1405 if(!committedPositions.count(node->sID) || (mapped-node->get_P()).norm()<=tolerance)
1406 {
1407 committedPositions[node->sID]=mapped;
1408 }
1409 }
1410 }
1411 }
1412 std::map<size_t,std::vector<NetPtr>> remappedOriginals;
1413 std::map<std::pair<size_t,size_t>,std::vector<NetPtr>> sharedBoundaryCorners;
1414 std::map<std::pair<std::vector<size_t>,size_t>,std::vector<NetPtr>> enclosedCorners;
1415 for(auto& plan : plans)
1416 {
1417 // Only fully shared old contour lineages may share an enclosed
1418 // replacement. Equal coordinates alone never establish ownership.
1419 std::vector<size_t> enclosureOwners;
1420 const bool enclosed=std::any_of(plan.contour.begin(),plan.contour.end(),[](const Record& record)
1421 {
1422 return std::get<6>(record).enclosureLayer.has_value();
1423 });
1424 if(enclosed)
1425 {
1426 for(size_t arc : plan.arcs)
1427 {
1428 enclosureOwners.push_back(root(arc));
1429 }
1430 // Canonicalize only cyclic start and traversal direction, not
1431 // the order of visits: different sheet connectivity must stay
1432 // distinct even if it uses the same set of old junction arcs.
1433 auto best=enclosureOwners;
1434 for(bool reverse : {false,true})
1435 {
1436 auto candidate=enclosureOwners;
1437 if(reverse)
1438 {
1439 std::reverse(candidate.begin(),candidate.end());
1440 }
1441 for(size_t start=0;start<candidate.size();++start)
1442 {
1443 best=std::min(best,candidate);
1444 std::rotate(candidate.begin(),candidate.begin()+1,candidate.end());
1445 }
1446 }
1447 enclosureOwners=std::move(best);
1448 }
1449 std::set<const LoopNodeType*> reused;
1450 for(const auto& record : plan.contour)
1451 {
1452 const auto& info=std::get<6>(record);
1453 const auto patch=plan.loop->periodicGlidePlane->getPatch(std::get<1>(record));
1454 const VectorDim position=plan.loop->periodicGlidePlane->referencePlane->globalPosition(std::get<0>(record));
1455 const VectorDim mapped=position+patch->shift;
1456 std::pair<std::shared_ptr<Edge>,std::shared_ptr<Edge>> edges;
1457 if(std::get<2>(record).first>=0)
1458 {
1459 edges.first=patch->edges().at(std::get<2>(record).first);
1460 if(std::get<2>(record).second>=0)
1461 {
1462 edges.second=patch->edges().at(std::get<2>(record).second);
1463 }
1464 }
1465 NodePtr chosen;
1466 NetPtr networkNode;
1467 if(info.originalVertex)
1468 {
1469 const auto old=plan.anchors.at(*info.originalVertex);
1470 networkNode=old->networkNode;
1471 if((committedPositions.at(networkNode->sID)-mapped).norm()>tolerance)
1472 {
1473 auto& alternatives=remappedOriginals[networkNode->sID];
1474 networkNode.reset();
1475 for(const auto& other : alternatives)
1476 {
1477 if((other->get_P()-mapped).norm()<=tolerance)
1478 {
1479 networkNode=other;
1480 break;
1481 }
1482 }
1483 if(!networkNode)
1484 {
1485 networkNode=this->networkNodes().create(mapped,old->networkNode->get_V(),old->networkNode->climbVelocityScalar,old->networkNode->velocityReduction());
1486 alternatives.push_back(networkNode);
1487 }
1488 }
1489 // Recompute confinement: sliding can change which edge
1490 // contains an anchor, especially at mixed corners.
1491 edges={nullptr,nullptr};
1492 for(const auto& edge : patch->edges())
1493 {
1494 if((edge->meshIntersection->snap(mapped)-mapped).norm()<=tolerance)
1495 {
1496 if(!edges.first)
1497 {
1498 edges.first=edge;
1499 }
1500 else
1501 {
1502 edges.second=edge;
1503 }
1504 }
1505 }
1506 const auto& previous=old->periodicPlaneEdge;
1507 if(old->networkNode==networkNode && old->periodicPlanePatch()==patch &&
1508 ((previous.first==edges.first && previous.second==edges.second) ||
1509 (previous.first==edges.second && previous.second==edges.first)))
1510 {
1511 chosen=old;
1512 }
1513 }
1514 else if(info.polygonLocation)
1515 {
1516 const auto& location=*info.polygonLocation;
1517 const size_t arcID=plan.arcs.at(location.segment.sourceIndex);
1518 const double t=arcs[arcID].reversed ? 1-location.parameter : location.parameter;
1519 auto& candidates=crossings[root(arcID)];
1520 for(const auto& candidate : candidates)
1521 {
1522 if(std::abs(candidate.parameter-t)<=1.e-8 && (candidate.node->get_P()-mapped).norm()<=tolerance)
1523 {
1524 networkNode=candidate.node;
1525 break;
1526 }
1527 }
1528 if(!networkNode)
1529 {
1530 networkNode=this->networkNodes().create(mapped);
1531 candidates.push_back({t,networkNode});
1532 }
1533 for(const auto& old : arcs[arcID].subdivisions)
1534 {
1535 if(old->networkNode==networkNode && old->periodicPlanePatch()==patch &&
1536 (old->get_P()-position).norm()<=tolerance && !danglingBoundaryLoopNodes.count(old.get()))
1537 {
1538 const auto& oldEdges=old->periodicPlaneEdge;
1539 if((oldEdges.first==edges.first && oldEdges.second==edges.second) ||
1540 (oldEdges.first==edges.second && oldEdges.second==edges.first))
1541 {
1542 chosen=old;
1543 break;
1544 }
1545 }
1546 }
1547 }
1548 if(!networkNode && info.enclosureLayer)
1549 {
1550 auto& candidates=enclosedCorners[{enclosureOwners,*info.enclosureLayer}];
1551 for(const auto& candidate : candidates)
1552 {
1553 if((candidate->get_P()-mapped).norm()<=tolerance)
1554 {
1555 networkNode=candidate;
1556 break;
1557 }
1558 }
1559 if(!networkNode)
1560 {
1561 networkNode=this->networkNodes().create(mapped);
1562 candidates.push_back(networkNode);
1563 }
1564 }
1565 // A boundary-following corner inherits a shared junction only
1566 // when BOTH ends of its boundary walk originate on shared arcs.
1567 // Matching coordinates alone must not join unrelated loops.
1568 if(!networkNode && info.onPhysicalBoundary && !plan.arcs.empty())
1569 {
1570 const size_t index=&record-plan.contour.data();
1571 std::optional<size_t> before,after;
1572 for(size_t step=1;step<plan.contour.size();++step)
1573 {
1574 if(!before)
1575 {
1576 const auto& metadata=std::get<6>(plan.contour[(index+plan.contour.size()-step)%plan.contour.size()]);
1577 if(metadata.polygonLocation)
1578 {
1579 before=root(plan.arcs.at(metadata.polygonLocation->segment.sourceIndex));
1580 }
1581 }
1582 if(!after)
1583 {
1584 const auto& metadata=std::get<6>(plan.contour[(index+step)%plan.contour.size()]);
1585 if(metadata.polygonLocation)
1586 {
1587 after=root(plan.arcs.at(metadata.polygonLocation->segment.sourceIndex));
1588 }
1589 }
1590 if(before && after)
1591 {
1592 break;
1593 }
1594 }
1595 if(before && after)
1596 {
1597 auto& candidates=sharedBoundaryCorners[std::minmax(*before,*after)];
1598 for(const auto& candidate : candidates)
1599 {
1600 if((candidate->get_P()-mapped).norm()<=tolerance)
1601 {
1602 networkNode=candidate;
1603 break;
1604 }
1605 }
1606 if(!networkNode)
1607 {
1608 networkNode=this->networkNodes().create(mapped);
1609 candidates.push_back(networkNode);
1610 }
1611 }
1612 }
1613 if(!networkNode)
1614 {
1615 networkNode=this->networkNodes().create(mapped);
1616 }
1617 if(!chosen || reused.count(chosen.get()))
1618 {
1619 chosen=this->loopNodes().create(plan.loop,networkNode,position,patch,edges);
1620 }
1621 if(info.originalVertex && (chosen->periodicPlaneEdge.first || plan.anchors.at(*info.originalVertex)->periodicContourVertex || promoted.count(plan.anchors.at(*info.originalVertex))))
1622 {
1623 promoted.insert(chosen);
1624 }
1625 reused.insert(chosen.get());
1626 plan.newNodes.push_back(chosen);
1627 }
1628 if(!plan.newNodes.empty() && plan.newNodes.size()<3)
1629 {
1630 throw std::runtime_error("Boundary update produced a degenerate retained contour.");
1631 }
1632 }
1633 // Prepare the edit journal before mutating any committed objects.
1634 using LinkMap=typename std::remove_reference_t<decltype(this->loopLinks())>;
1635 using Key=typename LoopLinkType::KeyType;
1636 std::vector<Key> obsolete,created;
1637 std::vector<std::tuple<NodePtr,NodePtr,std::shared_ptr<LoopType>>> additions;
1638 std::vector<typename LinkMap::node_type> removed;
1639 struct PositionState
1640 {
1641 NetPtr node;
1642 typename NetworkNodeType::SplineNodeType spline;
1643 const Simplex<dim,dim>* simplex;
1644 };
1645 std::vector<PositionState> savedPositions;
1646 for(const auto& entry : committedPositions)
1647 {
1648 const auto node=this->networkNodes().get(entry.first);
1649 savedPositions.push_back({node,static_cast<const typename NetworkNodeType::SplineNodeType&>(*node),node->p_Simplex});
1650 }
1651 for(const auto& plan : plans)
1652 {
1653 std::set<Key> desired,replacements;
1654 for(size_t i=0;i<plan.newNodes.size();++i)
1655 {
1656 const auto a=plan.newNodes[i];
1657 const auto b=plan.newNodes[(i+1)%plan.newNodes.size()];
1658 const auto key=LoopLinkType::loopLinkKey(plan.loop,a,b);
1659 if(!desired.insert(key).second)
1660 {
1661 throw std::runtime_error("Boundary contour repeats a loop-link identity.");
1662 }
1663 const auto found=this->loopLinks().find(key);
1664 if(found==this->loopLinks().end() || found->second.source!=a)
1665 {
1666 additions.emplace_back(a,b,plan.loop);
1667 replacements.insert(key);
1668 }
1669 }
1670 for(size_t i=0;i<plan.oldNodes.size();++i)
1671 {
1672 const auto a=plan.oldNodes[i];
1673 const auto b=plan.oldNodes[(i+1)%plan.oldNodes.size()];
1674 const auto key=LoopLinkType::loopLinkKey(plan.loop,a,b);
1675 if(!desired.count(key) || replacements.count(key))
1676 {
1677 obsolete.push_back(key);
1678 }
1679 }
1680 }
1681 struct SegmentState
1682 {
1683 std::shared_ptr<NetworkLinkType> link;
1684 VectorDim burgers;
1685 double norm;
1686 std::shared_ptr<SlipSystem> slip;
1687 };
1688 std::map<const NetworkLinkType*,SegmentState> segmentStates;
1689 const auto saveSegment=[&](const std::shared_ptr<NetworkLinkType>& link)
1690 {
1691 if(link)
1692 {
1693 segmentStates.try_emplace(link.get(),SegmentState{link,link->Burgers,link->BurgersNorm,link->_slipSystem});
1694 }
1695 };
1696 for(const auto& key : obsolete)
1697 {
1698 saveSegment(this->loopLinks().at(key).networkLink());
1699 }
1700 for(const auto& item : additions)
1701 {
1702 saveSegment(this->networkLinks().get(NetworkLinkType::getKey(std::get<0>(item)->networkNode,std::get<1>(item)->networkNode)));
1703 }
1704 removed.reserve(obsolete.size());
1705 created.reserve(additions.size());
1706 try
1707 {
1708 for(const auto& entry : committedPositions)
1709 {
1710 this->networkNodes().get(entry.first)->set_P(entry.second);
1711 if(checkpoint)
1712 {
1713 checkpoint();
1714 }
1715 }
1716 for(const auto& key : obsolete)
1717 {
1718 auto handle=this->loopLinks().extract(key);
1719 handle.mapped().detachTopology();
1720 removed.push_back(std::move(handle));
1721 if(checkpoint)
1722 {
1723 checkpoint();
1724 }
1725 }
1726 for(const auto& item : additions)
1727 {
1728 const auto& [a,b,loop]=item;
1729 this->connect(a,b,loop);
1730 created.push_back(LoopLinkType::loopLinkKey(loop,a,b));
1731 if(checkpoint)
1732 {
1733 checkpoint();
1734 }
1735 }
1736 }
1737 catch(...)
1738 {
1739 for(const auto& key : created)
1740 {
1741 this->loopLinks().erase(key);
1742 }
1743 for(auto& handle : removed)
1744 {
1745 handle.mapped().attachTopology();
1746 this->loopLinks().insert(std::move(handle));
1747 }
1748 for(const auto& entry : savedPositions)
1749 {
1750 static_cast<typename NetworkNodeType::SplineNodeType&>(*entry.node)=entry.spline;
1751 entry.node->p_Simplex=entry.simplex;
1752 }
1753 for(const auto& entry : savedPositions)
1754 {
1755 entry.node->updateGeometry();
1756 }
1757 for(const auto& entry : segmentStates)
1758 {
1759 const auto& state=entry.second;
1760 state.link->Burgers=state.burgers;
1761 state.link->BurgersNorm=state.norm;
1762 state.link->_slipSystem=state.slip;
1763 state.link->straight.updateGeometry();
1764 }
1765 throw;
1766 }
1767 for(const auto& node : promoted)
1768 {
1769 node->periodicContourVertex=true;
1770 }
1771 danglingBoundaryLoopNodes.clear();
1772 }
1773
1774 template <int dim>
1775 int DislocationNetwork<dim>::verboseDislocationNetwork=0;
1776
1777 template class DislocationNetwork<3>;
1778}
1779#endif
model::DefectiveCrystalConfigIO::loopLinks
const std::vector< DislocationLoopLinkIO< dim > > & loopLinks() const
Definition DefectiveCrystalConfigIO.cpp:151
model::DefectiveCrystalConfigIO::loops
const std::vector< DislocationLoopIO< dim > > & loops() const
Definition DefectiveCrystalConfigIO.cpp:124
model::DefectiveCrystalConfigIO::nodes
const std::vector< DislocationNodeIO< dim > > & nodes() const
Definition DefectiveCrystalConfigIO.cpp:100
model::DefectiveCrystalConfigIO::loopNodes
const std::vector< DislocationLoopNodeIO< dim > > & loopNodes() const
Definition DefectiveCrystalConfigIO.cpp:112
model::DislocationLoopNode::initFromFile
static void initFromFile(const std::string &)
Definition DislocationLoopNode.cpp:72
model::DislocationLoopNode::SplineNodeType
SplineNode< DislocationLoopNode< dim >, dim, 0, Hermite > SplineNodeType
Definition DislocationLoopNode.h:46
model::DislocationNetwork::VectorDim
TraitsType::VectorDim VectorDim
Definition DislocationNetwork.h:91
model::DislocationNetwork::networkRemesher
DislocationNetworkRemesh< dim > networkRemesher
Definition DislocationNetwork.h:126
model::DislocationNetwork::networkLengthPerSlipSystem
std::vector< std::tuple< double, double, double, double > > networkLengthPerSlipSystem() const
Definition DislocationNetwork.cpp:409
model::DislocationNetwork::updateBoundaryNodes
void updateBoundaryNodes(const std::function< void()> &checkpoint={})
Definition DislocationNetwork.cpp:1211
model::DislocationNetwork::subcyclingBins
const std::set< int > subcyclingBins
Definition DislocationNetwork.h:149
model::DislocationNetwork::executeSingleGlideStepDiscreteEvents
void executeSingleGlideStepDiscreteEvents(const long int &runID)
Definition DislocationNetwork.cpp:1082
model::DislocationNetwork::SimplexDim
MicrostructureBase< dim >::SimplexDim SimplexDim
Definition DislocationNetwork.h:97
model::DislocationNetwork::inclusions
const std::shared_ptr< InclusionMicrostructure< dim > > & inclusions() const
Definition DislocationNetwork.cpp:510
model::DislocationNetwork::VectorMSize
MicrostructureBase< dim >::VectorMSize VectorMSize
Definition DislocationNetwork.h:98
model::DislocationNetwork::displacement
VectorDim displacement(const VectorDim &, const NodeType *const, const ElementType *const, const SimplexDim *const) const override
Definition DislocationNetwork.cpp:523
model::DislocationNetwork::DislocationNetwork
DislocationNetwork(MicrostructureContainerType &mc)
Definition DislocationNetwork.cpp:25
model::DislocationNetwork::nodalVelocityConstraints
const Eigen::Matrix< double, Eigen::Dynamic, dim > nodalVelocityConstraints
Definition DislocationNetwork.h:148
model::DislocationNetwork::outputDislocationDensityPerSlipSystem
const bool outputDislocationDensityPerSlipSystem
Definition DislocationNetwork.h:142
model::DislocationNetwork::junctionsMaker
DislocationJunctionFormation< dim > junctionsMaker
Definition DislocationNetwork.h:127
model::DislocationNetwork::slipSystemAveragePlasticDistortion
std::map< std::pair< int, int >, double > slipSystemAveragePlasticDistortion() const
Definition DislocationNetwork.cpp:356
model::DislocationNetwork::NodeType
MicrostructureBase< dim >::NodeType NodeType
Definition DislocationNetwork.h:96
model::DislocationNetwork::VectorISize
MicrostructureBase< dim >::VectorISize VectorISize
Definition DislocationNetwork.h:99
model::DislocationNetwork::networkLength
std::tuple< double, double, double, double > networkLength() const
Definition DislocationNetwork.cpp:462
model::DislocationNetwork::NdofXnode
static constexpr int NdofXnode
Definition DislocationNetwork.h:102
model::DislocationNetwork::glideSolver
std::shared_ptr< DislocationGlideSolverBase< dim > > glideSolver
Definition DislocationNetwork.h:115
model::DislocationNetwork::ddBase
DefectiveCrystalBase< dim > & ddBase
Definition DislocationNetwork.h:125
model::DislocationNetwork::averagePlasticDistortion
MatrixDim averagePlasticDistortion() const override
Definition DislocationNetwork.cpp:337
model::DislocationNetwork::output
void output(DefectiveCrystalConfigIO< dim > &configIO, DefectiveCrystalAuxIO< dim > &auxIO, std::ofstream &f_file, std::ofstream &F_labels) const override
Definition DislocationNetwork.cpp:1104
model::DislocationNetwork::stochasticForceGenerator
std::shared_ptr< StochasticForceGenerator > stochasticForceGenerator
Definition DislocationNetwork.h:132
model::DislocationNetwork::contract
bool contract(std::shared_ptr< NetworkNodeType > nA, std::shared_ptr< NetworkNodeType > nB)
Definition DislocationNetwork.cpp:516
model::DislocationNetwork::setConfiguration
void setConfiguration(const DefectiveCrystalConfigIO< dim > &)
Definition DislocationNetwork.cpp:231
model::DislocationNetwork::stress
MatrixDim stress(const VectorDim &, const NodeType *const, const ElementType *const ele, const SimplexDim *const guess) const override
Definition DislocationNetwork.cpp:634
model::DislocationNetwork::averagePlasticStrain
MatrixDim averagePlasticStrain() const
Definition DislocationNetwork.cpp:348
model::DislocationNetwork::mobileConcentration
VectorMSize mobileConcentration(const VectorDim &, const NodeType *const, const ElementType *const, const SimplexDim *const) const override
Definition DislocationNetwork.cpp:564
model::DislocationNetwork::MatrixDim
TraitsType::MatrixDim MatrixDim
Definition DislocationNetwork.h:93
model::DislocationNetwork::averagePlasticStrainRate
MatrixDim averagePlasticStrainRate() const
Definition DislocationNetwork.cpp:401
model::DislocationNetwork::timeStepper
DDtimeStepper< dim > timeStepper
Definition DislocationNetwork.h:131
model::DislocationNetwork::immobileClusters
VectorISize immobileClusters(const VectorDim &, const NodeType *const, const ElementType *const, const SimplexDim *const) const override
Definition DislocationNetwork.cpp:627
model::DislocationNetwork::moveNodes
void moveNodes(const double &dt_in)
Definition DislocationNetwork.cpp:1006
model::DislocationNetwork::inelasticDisplacementRate
VectorDim inelasticDisplacementRate(const VectorDim &, const NodeType *const, const ElementType *const, const SimplexDim *const) const override
Definition DislocationNetwork.cpp:656
model::DislocationNetwork::MicrostructureContainerType
MicrostructureContainer< dim > MicrostructureContainerType
Definition DislocationNetwork.h:94
model::DislocationNetwork::storeSingleGlideStepDiscreteEvents
void storeSingleGlideStepDiscreteEvents(const long int &runID)
Definition DislocationNetwork.cpp:1072
model::DislocationNetwork::outputPlasticDistortionPerSlipSystem
const bool outputPlasticDistortionPerSlipSystem
Definition DislocationNetwork.h:141
model::DislocationNetwork::nodeContractor
DislocationNodeContraction< dim > nodeContractor
Definition DislocationNetwork.h:130
model::DislocationNetwork::averagePlasticDistortionRate
MatrixDim averagePlasticDistortionRate() const override
Definition DislocationNetwork.cpp:390
model::DislocationNetwork::addConfiguration
void addConfiguration(const DefectiveCrystalConfigIO< dim > &)
Definition DislocationNetwork.cpp:88
model::DislocationNetwork::danglingBoundaryLoopNodes
std::set< const LoopNodeType * > danglingBoundaryLoopNodes
Definition DislocationNetwork.h:145
model::DislocationNetwork::crossSlipMaker
DislocationCrossSlip< dim > crossSlipMaker
Definition DislocationNetwork.h:129
model::DislocationNetwork::climbSolver
std::shared_ptr< DislocationClimbSolverBase< dim > > climbSolver
Definition DislocationNetwork.h:116
model::DislocationNetwork::getDt
double getDt() const override
Definition DislocationNetwork.cpp:680
model::DislocationNetwork::ElementType
MicrostructureBase< dim >::ElementType ElementType
Definition DislocationNetwork.h:95
model::DislocationNetwork::averageStress
MatrixDim averageStress() const override
Definition DislocationNetwork.cpp:554
model::DislocationNetwork::initializeConfiguration
void initializeConfiguration(const DefectiveCrystalConfigIO< dim > &configIO, const std::ofstream &f_file, const std::ofstream &F_labels) override
Definition DislocationNetwork.cpp:70
model::DislocationNetwork::crossSlipModel
const std::shared_ptr< BaseCrossSlipModel< dim > > crossSlipModel
Definition DislocationNetwork.h:128
model::DislocationSegment::initFromFile
static void initFromFile(const std::string &)
Definition DislocationSegment.cpp:234
model::LoopNetwork< DislocationNetwork< dim > >::loops
const LoopContainerType & loops() const
Definition LoopNetwork.h:98
model::LoopNetwork< DislocationNetwork< dim > >::connect
void connect(const SharedLoopNodePtrType &n0, const SharedLoopNodePtrType &n1, const std::shared_ptr< LoopType > &tempLoop)
Definition LoopNetwork.h:211
model::LoopNetwork< DislocationNetwork< dim > >::networkNodes
const NetworkNodeContainerType & networkNodes() const
Definition LoopNetwork.h:118
model::LoopNetwork< DislocationNetwork< dim > >::loopNodes
const LoopNodeContainerType & loopNodes() const
Definition LoopNetwork.h:108
model::LoopNetwork< DislocationNetwork< dim > >::networkLinks
NetworkLinkContainerType & networkLinks()
Definition LoopNetwork.h:139
model::LoopNetwork< DislocationNetwork< dim > >::insertLoop
void insertLoop(const std::shared_ptr< LoopType > &loop, const std::vector< std::shared_ptr< LoopNodeType > > nodes)
Definition LoopNetwork.h:246
model::StaticID< DislocationLoopNode< dim > >::set_count
static void set_count(const size_t &newCount)
Definition StaticID.h:65
model::StaticID< DislocationNode< dim > >::force_count
static void force_count(const size_t &newCount)
Definition StaticID.h:79
model::TextFileParser
Definition TextFileParser.h:123
model::TextFileParser::readString
std::string readString(const std::string &key, const bool &verbose=false)
Definition TextFileParser.h:202
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
model::count
size_t InstanceCount< Derived >::count
Definition InstanceCount.h:70
model::retainPeriodicContourBends
void retainPeriodicContourBends(const Nodes &nodes, const double tolerance)
Definition PeriodicContourTopology.h:13
model::magentaColor
static std::string magentaColor
Definition TerminalColors.h:31
model::retracesContourSegment
bool retracesContourSegment(const Nodes &nodes, const std::vector< size_t > &anchors, const double tolerance)
Definition PeriodicContourTopology.h:46
model::PolygonPatchIntersectionRecord
std::tuple< Eigen::Matrix< double, dim-1, 1 >, Eigen::Matrix< double, dim, 1 >, std::pair< short int, short int >, std::map< Eigen::Matrix< double, dim, 1 >, std::set< std::pair< short int, short int > >, CompareVectorsByComponent< double, dim, float > >, size_t, const T *const, PolygonPatchMetadata< T > > PolygonPatchIntersectionRecord
Definition PolygonPatchIntersection.h:68
model::defaultColor
static std::string defaultColor
Definition TerminalColors.h:17
model::kd
bool kd(const int &i, const int &j)
Definition ElasticDeformation.cpp:21
model::GlidePlaneKey
LatticePlaneKey< dim > GlidePlaneKey
Definition GlidePlaneKey.h:20
model::DefectiveCrystalAuxIO::quadraturePoints
const std::vector< DislocationQuadraturePointIO< dim > > & quadraturePoints() const
Definition DefectiveCrystalAuxIO.cpp:77
model::DislocationClimbSolverFactory::getClimbSolver
static std::shared_ptr< DislocationClimbSolverBase< dim > > getClimbSolver(const DislocationNetwork< dim > &DN, const std::string &solverType)
Definition DislocationClimbSolverFactory.cpp:64
model::DislocationFieldBase::initFromFile
static void initFromFile(const std::string &fileName)
dim x dim identity matrix
Definition DislocationFieldBase.cpp:20
model::DislocationGlideSolverFactory::getGlideSolver
static std::shared_ptr< DislocationGlideSolverBase< dim > > getGlideSolver(const DislocationNetwork< dim > &DN, const std::string &solverType)
Definition DislocationGlideSolverFactory.cpp:31
model::DislocationNucleation::bulkNucleate
void bulkNucleate(const int &bulkNucleationModel)
Definition DislocationNucleation.cpp:35
model::DislocationNucleation::surfaceNucleate
void surfaceNucleate(const int &surfaceNucleationModel)
Definition DislocationNucleation.cpp:76
model::MicrostructureBase::pointGrains
std::set< const Grain< dim > * > pointGrains(const VectorDim &x, const NodeType *const node, const ElementType *const ele, const SimplexDim *const guess) const
Definition src/DefectiveCrystal/MicrostructureBase.cpp:71
model::MicrostructureBase::MicrostructureBase
MicrostructureBase(const std::string &tag_in, MicrostructureContainer< dim > &m_in)
Definition src/DefectiveCrystal/MicrostructureBase.cpp:25
model::MicrostructureBase::microstructures
MicrostructureContainerType & microstructures
Definition MicrostructureBase.h:47
parameter
void parameter(std::string &text, const std::string &key, const std::string &value)
Definition testClusterPeriodic.cpp:23
end
int end[2]
Definition triangle.cpp:11159
edges
m edges
Definition triangle.cpp:15849
j
int j
Definition triangle.cpp:3300
t0
REAL t0
Definition triangle.cpp:5171
intersect
enum locateresult intersect
Definition triangle.cpp:8237
index
int index
Definition triangle.cpp:12439
axis
int axis
Definition triangle.cpp:9283
k
int k
Definition triangle.cpp:3300
info
void info()
Definition triangle.cpp:1548
i
int i
Definition triangle.cpp:3300
key
newbad key
Definition triangle.cpp:7029
temp
vertex temp
Definition triangle.cpp:9214
poly
b poly
Definition triangle.cpp:3303
first
vertex first
Definition triangle.cpp:8231
Generated by doxygen 1.15.0