← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/atc/DislocationExtractor.h

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2013-11-20 22:41:36 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20131120224136-tzx7leh606fqnckm
Tags: 0~20131119.git7162cf0-1
* [e65b919] Imported Upstream version 0~20131119.git7162cf0
* [f7bddd4] Fix some problems, introduced by upstream recently.
* [3616dfc] Use wrap-and-sort script.
* [7e92030] Ignore quilt dir

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/atc/DislocationExtractor.h
 
1
///////////////////////////////////////////////////////////////////////////////
 
2
//
 
3
//  Copyright 2010, Alexander Stukowski
 
4
//  All rights reserved. See README.txt for more information.
 
5
//
 
6
///////////////////////////////////////////////////////////////////////////////
 
7
 
 
8
#ifndef __DXA_DISLOCATION_EXTRACTOR_H
 
9
#define __DXA_DISLOCATION_EXTRACTOR_H
 
10
 
 
11
// Include the LAMMPS headers that we need.
 
12
#include "atom.h"
 
13
#include "comm.h"
 
14
#include "error.h"
 
15
#include "domain.h"
 
16
#include "update.h"
 
17
#include "neigh_list.h"
 
18
using namespace LAMMPS_NS;
 
19
 
 
20
// This is must be defined when we build the LAMMPS module (and not the DXA standalone tool).
 
21
#define DXA_LAMMPS_MODULE
 
22
 
 
23
// Enables additional sanity checks.
 
24
#define DEBUG_DISLOCATIONS
 
25
 
 
26
// In LAMMPS, we always use double precision floating-point numbers.
 
27
#define USE_DOUBLE_PRECISION_FP
 
28
 
 
29
#include "../src/dxa/DXATracing.h"
 
30
#include "../src/util/FILEStream.h"
 
31
 
 
32
/**
 
33
 *
 
34
 */
 
35
class DXADislocationExtractor : private DXATracing, Pointers
 
36
{
 
37
public:
 
38
 
 
39
  /// Constructor.
 
40
  ///
 
41
  /// Expects a pointer to the global LAMMPS object.
 
42
  ///
 
43
  /// Exact mode flag:
 
44
  ///
 
45
  /// Enables serial processing, i.e. the master processor performs the complete analysis.
 
46
  /// This allows for proper handling of periodic boundary conditions.
 
47
  DXADislocationExtractor(LAMMPS* lmp, bool _exactMode = false) : DXATracing(msgLoggerStream, verboseLoggerStream), Pointers(lmp)
 
48
  {
 
49
    exactMode = _exactMode;
 
50
 
 
51
    // Connect internal logging streams to global LAMMPS streams.
 
52
    if(comm->me == 0) {
 
53
      msgLoggerStream.rdbuf()->setFiles(screen, logfile);
 
54
      verboseLoggerStream.rdbuf()->setFiles(screen, logfile);
 
55
    }
 
56
    this->processor = comm->me;
 
57
  }
 
58
 
 
59
  /// Extracts all dislocations from the current simulation state.
 
60
  ///
 
61
  /// First copies atoms from LAMMPS array to internal memory.
 
62
  /// Also adopts the nearest neighbor lists from LAMMPS.
 
63
  ///
 
64
  /// Parameters:
 
65
  ///
 
66
  ///      neighborList   A full neighbor list provided by LAMMPS, which contains neighbors of ghost atoms.
 
67
  ///      cnaCutoff      Cutoff radius to be used for the common neighbor analysis.
 
68
  ///
 
69
  void extractDislocations(NeighList* neighborList, double cnaCutoff, int lineSmoothingLevel = DEFAULT_LINE_SMOOTHING_LEVEL, int lineCoarseningLevel = DEFAULT_LINE_COARSENING_LEVEL)
 
70
  {
 
71
    // Check parameters.
 
72
    if(cnaCutoff <= 0.0) lammps->error->all(FLERR,"Common neighbor analysis cutoff radius must be positive.");
 
73
 
 
74
    // Release any memory allocated during last call.
 
75
    cleanup();
 
76
 
 
77
    // Setup simulation cell.
 
78
    if(exactMode) {
 
79
      pbc[0] = domain->periodicity[0] != 0;
 
80
      pbc[1] = domain->periodicity[1] != 0;
 
81
      pbc[2] = domain->periodicity[2] != 0;
 
82
    }
 
83
    else {
 
84
      pbc[0] = pbc[1] = pbc[2] = false;
 
85
    }
 
86
    simulationCell(0,0) = domain->h[0];
 
87
    simulationCell(1,1) = domain->h[1];
 
88
    simulationCell(2,2) = domain->h[2];
 
89
    simulationCell(1,2) = domain->h[3];
 
90
    simulationCell(0,2) = domain->h[4];
 
91
    simulationCell(0,1) = domain->h[5];
 
92
    simulationCell(1,0) = simulationCell(2,0) = simulationCell(2,1) = 0;
 
93
    simulationCellOrigin.X = domain->boxlo[0];
 
94
    simulationCellOrigin.Y = domain->boxlo[1];
 
95
    simulationCellOrigin.Z = domain->boxlo[2];
 
96
    setCNACutoff(cnaCutoff);
 
97
    setupSimulationCell(cnaCutoff);
 
98
    timestep = update->ntimestep;
 
99
 
 
100
    // Transfer LAMMPS atoms to internal array.
 
101
    if(exactMode)
 
102
      transferAllLAMMPSAtoms();
 
103
    else
 
104
      transferLocalLAMMPSAtoms(neighborList);
 
105
 
 
106
    if(exactMode == false || processor == 0) {
 
107
 
 
108
      // Perform common neighbor analysis to identify crystalline atoms.
 
109
      performCNA();
 
110
 
 
111
      // Order the neighbor lists of crystalline atoms.
 
112
      orderCrystallineAtoms();
 
113
 
 
114
      // Cluster crystalline atoms.
 
115
      clusterAtoms();
 
116
 
 
117
#ifdef PROCESSED_ATOMS
 
118
      // Enable this code block to get the processed atoms.
 
119
      {
 
120
        ostringstream ss;
 
121
        ss << "proc_" << comm->me << ".atoms.dump";
 
122
        ofstream stream(ss.str().c_str());
 
123
        writeAtomsDumpFile(stream);
 
124
      }
 
125
#endif
 
126
 
 
127
      // Create the nodes of the interface mesh.
 
128
      createInterfaceMeshNodes();
 
129
 
 
130
      // Create the facets of the interface mesh.
 
131
      createInterfaceMeshFacets();
 
132
 
 
133
#ifdef DEBUG_DISLOCATIONS
 
134
      // Check the generated mesh.
 
135
      validateInterfaceMesh();
 
136
#endif
 
137
 
 
138
      // Generate and advance Burgers circuits on the interface mesh.
 
139
      traceDislocationSegments();
 
140
 
 
141
      // Smooth dislocation lines.
 
142
      smoothDislocationSegments(lineSmoothingLevel, lineCoarseningLevel);
 
143
    }
 
144
 
 
145
    if(exactMode) {
 
146
      // Wrap segments at simulation cell boundaries.
 
147
      wrapDislocationSegments();
 
148
      // Distribute extracted dislocation segments to back all processors.
 
149
      broadcastDislocations();
 
150
    }
 
151
 
 
152
    // Clip dislocation lines at processor domain.
 
153
    Point3 subOrigin;
 
154
    Matrix3 subCell;
 
155
        if(domain->triclinic) {
 
156
      subOrigin = simulationCellOrigin + simulationCell * Vector3(domain->sublo_lamda);
 
157
      subCell.setColumn(0, simulationCell.column(0) * (domain->subhi_lamda[0] - domain->sublo_lamda[0]));
 
158
      subCell.setColumn(1, simulationCell.column(1) * (domain->subhi_lamda[1] - domain->sublo_lamda[1]));
 
159
      subCell.setColumn(2, simulationCell.column(2) * (domain->subhi_lamda[2] - domain->sublo_lamda[2]));
 
160
        } else {
 
161
      subOrigin = Point3(domain->sublo);
 
162
      subCell.setColumn(0, Vector3(domain->subhi[0] - domain->sublo[0], 0, 0));
 
163
      subCell.setColumn(1, Vector3(0, domain->subhi[1] - domain->sublo[1], 0));
 
164
      subCell.setColumn(2, Vector3(0, 0, domain->subhi[2] - domain->sublo[2]));
 
165
        }
 
166
    clipDislocationLines(subOrigin, subCell);
 
167
 
 
168
#ifdef OUTPUT_SEGMENTS
 
169
    // Enable this code block to see the extracted dislocation segments on each processor.
 
170
    {
 
171
      ostringstream ss;
 
172
      ss << "proc_" << comm->me << ".dislocations.vtk";
 
173
      ofstream stream(ss.str().c_str());
 
174
      writeDislocationsVTKFile(stream);
 
175
    }
 
176
#endif
 
177
  }
 
178
 
 
179
  /// Returns the list of extracted dislocation segments.
 
180
  const vector<DislocationSegment*>& getSegments() const { return this->segments; }
 
181
 
 
182
protected:
 
183
 
 
184
  /// Copies local atoms and neighbor lists from LAMMPS to internal array.
 
185
  void transferLocalLAMMPSAtoms(NeighList* neighborList)
 
186
  {
 
187
    // Allocate internal atoms array.
 
188
    int numTotalAtoms = atom->nlocal + atom->nghost;
 
189
    inputAtoms.resize(numTotalAtoms);
 
190
    numLocalInputAtoms = numTotalAtoms;
 
191
 
 
192
    // Transfer local and ghost atoms to internal array.
 
193
    for(int i = 0; i < numTotalAtoms; i++) {
 
194
      InputAtom& a = inputAtoms[i];
 
195
      a.tag = i + 1;
 
196
      a.pos.X = atom->x[i][0];
 
197
      a.pos.Y = atom->x[i][1];
 
198
      a.pos.Z = atom->x[i][2];
 
199
      a.flags = 0;
 
200
      a.cluster = NULL;
 
201
      a.numNeighbors = 0;
 
202
      a.setFlag(ATOM_IS_LOCAL_ATOM);
 
203
    }
 
204
 
 
205
    // Adopt nearest neighbor lists from LAMMPS.
 
206
    if(neighborList->ghostflag) {
 
207
      DISLOCATIONS_ASSERT(neighborList->inum + neighborList->gnum == numTotalAtoms);
 
208
      int tt = 0;
 
209
      for(int ii = 0; ii < numTotalAtoms; ii++) {
 
210
        int i = neighborList->ilist[ii];
 
211
        DISLOCATIONS_ASSERT(i < numTotalAtoms);
 
212
        InputAtom& ai = inputAtoms[i];
 
213
        int* jlist = neighborList->firstneigh[i];
 
214
        int jnum = neighborList->numneigh[i];
 
215
        if(ii >= neighborList->inum)
 
216
          tt += jnum;
 
217
        for(int jj = 0; jj < jnum; jj++) {
 
218
          int j = jlist[jj];
 
219
          j &= NEIGHMASK;
 
220
          DISLOCATIONS_ASSERT(j >= 0 && j < numTotalAtoms);
 
221
          if(j < i) continue;
 
222
          InputAtom& aj = inputAtoms[j];
 
223
          if(DistanceSquared(aj.pos, ai.pos) >= cnaCutoff * cnaCutoff) continue;
 
224
          if(ai.numNeighbors == MAX_ATOM_NEIGHBORS)
 
225
            raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", ai.tag, MAX_ATOM_NEIGHBORS);
 
226
          ai.addNeighbor(&aj);
 
227
          if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
 
228
            raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
 
229
          aj.addNeighbor(&ai);
 
230
        }
 
231
      }
 
232
    }
 
233
    else {
 
234
      for(int ii = 0; ii < neighborList->inum; ii++) {
 
235
        int i = neighborList->ilist[ii];
 
236
        InputAtom& ai = inputAtoms[i];
 
237
        int* jlist = neighborList->firstneigh[i];
 
238
        int jnum = neighborList->numneigh[i];
 
239
        for(int jj = 0; jj < jnum; jj++) {
 
240
          int j = jlist[jj];
 
241
          j &= NEIGHMASK;
 
242
          DISLOCATIONS_ASSERT(j >= 0 && j < numTotalAtoms);
 
243
          if(j < i) continue;
 
244
          InputAtom& aj = inputAtoms[j];
 
245
          if(DistanceSquared(aj.pos, ai.pos) < cnaCutoff * cnaCutoff) {
 
246
            if(ai.numNeighbors == MAX_ATOM_NEIGHBORS)
 
247
              raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", ai.tag, MAX_ATOM_NEIGHBORS);
 
248
            ai.addNeighbor(&aj);
 
249
            if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
 
250
              raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
 
251
            aj.addNeighbor(&ai);
 
252
          }
 
253
          if(j >= atom->nlocal) {
 
254
            for(int kk = 0; kk < jnum; kk++) {
 
255
              int k = jlist[kk];
 
256
              k &= NEIGHMASK;
 
257
              DISLOCATIONS_ASSERT(k < (int)inputAtoms.size());
 
258
              if(k == j) continue;
 
259
              if(k < atom->nlocal) continue;
 
260
              InputAtom& ak = inputAtoms[k];
 
261
              if(DistanceSquared(aj.pos, ak.pos) >= cnaCutoff * cnaCutoff) continue;
 
262
              if(aj.hasNeighbor(&ak)) continue;
 
263
              if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
 
264
                raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
 
265
              aj.addNeighbor(&ak);
 
266
            }
 
267
          }
 
268
        }
 
269
      }
 
270
    }
 
271
  }
 
272
 
 
273
  /// Copies atoms from LAMMPS to internal array and gathers all atoms on the master processor.
 
274
  void transferAllLAMMPSAtoms()
 
275
  {
 
276
    // Determine maximum number of local atoms a single processor has.
 
277
    int nlocalatoms_max;
 
278
    MPI_Reduce(&atom->nlocal, &nlocalatoms_max, 1, MPI_INT, MPI_MAX, 0, world);
 
279
 
 
280
      if(processor == 0) {
 
281
      // Allocate internal atoms array.
 
282
      inputAtoms.resize(atom->natoms);
 
283
      numLocalInputAtoms = atom->natoms;
 
284
      vector<InputAtom>::iterator currentAtom = inputAtoms.begin();
 
285
 
 
286
      // Copy local atoms of master processor into internal array.
 
287
      for(int i = 0; i < atom->nlocal; i++, ++currentAtom) {
 
288
        currentAtom->tag = currentAtom - inputAtoms.begin() + 1;
 
289
        currentAtom->pos.X = atom->x[i][0];
 
290
        currentAtom->pos.Y = atom->x[i][1];
 
291
        currentAtom->pos.Z = atom->x[i][2];
 
292
        currentAtom->flags = 0;
 
293
        currentAtom->cluster = NULL;
 
294
        currentAtom->numNeighbors = 0;
 
295
        currentAtom->setFlag(ATOM_IS_LOCAL_ATOM);
 
296
      }
 
297
 
 
298
      if(comm->nprocs > 1) {
 
299
        // Allocate receive buffer.
 
300
        vector<double> buffer(nlocalatoms_max * 3);
 
301
 
 
302
        // Receive atoms from other processors.
 
303
        for(int iproc = 1; iproc < comm->nprocs; iproc++) {
 
304
          MPI_Status status;
 
305
          MPI_Recv(buffer.empty() ? NULL : &buffer.front(), nlocalatoms_max * 3, MPI_DOUBLE, iproc, 0, world, &status);
 
306
          int ndoubles;
 
307
          MPI_Get_count(&status, MPI_DOUBLE, &ndoubles);
 
308
          int nReceived = ndoubles / 3;
 
309
 
 
310
          vector<double>::const_iterator data = buffer.begin();
 
311
          for(int i = 0; i < nReceived; i++, ++currentAtom) {
 
312
            currentAtom->tag = currentAtom - inputAtoms.begin() + 1;
 
313
            currentAtom->pos.X = *data++;
 
314
            currentAtom->pos.Y = *data++;
 
315
            currentAtom->pos.Z = *data++;
 
316
            currentAtom->flags = 0;
 
317
            currentAtom->cluster = NULL;
 
318
            currentAtom->numNeighbors = 0;
 
319
            currentAtom->setFlag(ATOM_IS_LOCAL_ATOM);
 
320
          }
 
321
        }
 
322
      }
 
323
        DISLOCATIONS_ASSERT(currentAtom == inputAtoms.end());
 
324
    }
 
325
    else {
 
326
      // Allocate and fill send buffer.
 
327
      vector<double> buffer(atom->nlocal * 3);
 
328
      vector<double>::iterator data = buffer.begin();
 
329
      for(int i = 0; i < atom->nlocal; i++) {
 
330
        *data++ = atom->x[i][0];
 
331
        *data++ = atom->x[i][1];
 
332
        *data++ = atom->x[i][2];
 
333
      }
 
334
      // Send local atom coordinates to master proc.
 
335
      MPI_Send(buffer.empty() ? NULL : &buffer.front(), buffer.size(), MPI_DOUBLE, 0, 0, world);
 
336
    }
 
337
 
 
338
      // Make sure all input atoms are wrapped at periodic boundary conditions.
 
339
      wrapInputAtoms(NULL_VECTOR);
 
340
 
 
341
      // Build nearest neighbor lists.
 
342
      buildNearestNeighborLists();
 
343
  }
 
344
 
 
345
  /// This structure is used to communicate dislocation segments between processors.
 
346
  struct DislocationSegmentComm {
 
347
    LatticeVector burgersVector;
 
348
    Vector3 burgersVectorWorld;
 
349
    int numPoints;
 
350
  };
 
351
 
 
352
  /// After the master processor has extracted all dislocation segments, this function broadcasts the dislocations
 
353
  /// back to all other processor.
 
354
  void broadcastDislocations()
 
355
  {
 
356
    if(comm->nprocs == 1) return;  // Nothing to do in serial mode.
 
357
 
 
358
    // Broadcast number of segments.
 
359
    int numSegments = segments.size();
 
360
    MPI_Bcast(&numSegments, 1, MPI_INT, 0, world);
 
361
 
 
362
    // Allocate send/receive buffer for dislocation segments.
 
363
    vector<DislocationSegmentComm> segmentBuffer(numSegments);
 
364
 
 
365
    // The total number of line points (sum of all segments).
 
366
    int numLinePoints = 0;
 
367
 
 
368
    if(processor == 0) {
 
369
      // Fill segment send buffer.
 
370
      vector<DislocationSegmentComm>::iterator sendItem = segmentBuffer.begin();
 
371
      for(vector<DislocationSegment*>::const_iterator segment = segments.begin(); segment != segments.end(); ++segment, ++sendItem) {
 
372
        sendItem->burgersVector = (*segment)->burgersVector;
 
373
        sendItem->burgersVectorWorld = (*segment)->burgersVectorWorld;
 
374
        sendItem->numPoints = (*segment)->line.size();
 
375
        numLinePoints += (*segment)->line.size();
 
376
      }
 
377
    }
 
378
    // Broadcast segments.
 
379
    MPI_Bcast(segmentBuffer.empty() ? NULL : &segmentBuffer.front(), segmentBuffer.size() * sizeof(segmentBuffer[0]), MPI_CHAR, 0, world);
 
380
 
 
381
    if(processor != 0) {
 
382
      // Extract segments from receive buffer.
 
383
      segments.reserve(segmentBuffer.size());
 
384
      for(vector<DislocationSegmentComm>::const_iterator receiveItem = segmentBuffer.begin(); receiveItem != segmentBuffer.end(); ++receiveItem) {
 
385
        DislocationSegment* newSegment = segmentPool.construct(receiveItem->burgersVector, receiveItem->burgersVectorWorld);
 
386
        newSegment->index = segments.size() + 1;
 
387
        segments.push_back(newSegment);
 
388
        numLinePoints += receiveItem->numPoints;
 
389
      }
 
390
    }
 
391
 
 
392
    // Allocate send/receive buffer for dislocation points.
 
393
    vector<Point3> pointBuffer(numLinePoints);
 
394
 
 
395
    if(processor == 0) {
 
396
      // Fill point send buffer.
 
397
      vector<Point3>::iterator sendItem = pointBuffer.begin();
 
398
      for(vector<DislocationSegment*>::const_iterator segment = segments.begin(); segment != segments.end(); ++segment) {
 
399
        for(deque<Point3>::const_iterator p = (*segment)->line.begin(); p != (*segment)->line.end(); ++p, ++sendItem)
 
400
          *sendItem = *p;
 
401
      }
 
402
      DISLOCATIONS_ASSERT(sendItem == pointBuffer.end());
 
403
    }
 
404
    // Broadcast segments.
 
405
    MPI_Bcast(pointBuffer.empty() ? NULL : &pointBuffer.front(), pointBuffer.size() * sizeof(pointBuffer[0]), MPI_CHAR, 0, world);
 
406
 
 
407
    if(processor != 0) {
 
408
      // Extract points from receive buffer.
 
409
      vector<DislocationSegment*>::const_iterator segment = segments.begin();
 
410
      vector<Point3>::const_iterator p = pointBuffer.begin();
 
411
      for(vector<DislocationSegmentComm>::const_iterator receiveItem = segmentBuffer.begin(); receiveItem != segmentBuffer.end(); ++receiveItem, ++segment) {
 
412
        (*segment)->line.assign(p, p + receiveItem->numPoints);
 
413
        p += receiveItem->numPoints;
 
414
      }
 
415
      DISLOCATIONS_ASSERT(p == pointBuffer.end());
 
416
    }
 
417
  }
 
418
 
 
419
protected:
 
420
 
 
421
  /// Pointer to the main LAMMPS object.
 
422
  LAMMPS* lammps;
 
423
 
 
424
  /// Enables serial processing, which provides exact handling of periodic boundary conditions.
 
425
  bool exactMode;
 
426
 
 
427
  /// The output stream to which log messages are sent.
 
428
  stdio_osyncstream msgLoggerStream;
 
429
 
 
430
  /// The output stream to which verbose log messages are sent.
 
431
  stdio_osyncstream verboseLoggerStream;;
 
432
};
 
433
 
 
434
#endif // __DXA_DISLOCATION_EXTRACTOR_H