1. add a non-singular displacement field method to generate any (7c451d47) · Commits · DeAn Wei / prodata

include/Home.h

+12 −0

Original line number	Diff line number	Diff line
		@@ -99,6 +99,11 @@ typedef enum{
		FTYPE_MAX
		}FTYPE_t;

		typedef enum{
		INT_DATA = 0,
		DOUBLE_DATA,
		STRING_DATA
		}DATA_TYPE;
		/*
		* Define a structure to hold a command line option's id (type),
		* name, the shortest possible unique abbreviation of the option
		@@ -111,8 +116,15 @@ typedef struct {
		int optPaired;
		} Option_t;

		typedef struct {
		int type;
		string name;
		string val;
		}Variable_t;

		typedef struct {
		vector<string> variables;
		vector<Variable_t> auxData;
		vector<vector<double> > data;
		} Table_t;

include/ProDataIO.h

+1 −0

Original line number	Diff line number	Diff line
		@@ -31,5 +31,6 @@ int ReadDumpFile(const string &file, Dump_t &dum);
		int ReadDataFromMDLogFile(const vector<string> &files, LineList_t &list);
		int WriteDumpFile(const string &file, Dump_t &dum, int precision = 6);
		int MGToLMPDataFile(const string &file, Dump_t &dum, int precision = 10);
		bool ReadLMPFile(const string file, Dump_t &dum);

		#endif

src/GenerateDislocation.cpp

+166 −7

Original line number	Diff line number	Diff line
		@@ -16,7 +16,7 @@ void GenerateScrewDislocation(InArgs_t *inArgs)
		Dump_t dum;
		string posName("pos"), lineName("line"), gDirName("gdir"), burgMagName("burgMag");
		string screwTypeName("stype"), multiDisName("multi");
		int screwType = 0, nPairs = 1;
		int screwType = 0, nPairs = 0;
		double rsize = 0, delta;

		real8 pos[3] = {0.0, 0.0, 0.0};
		@@ -73,7 +73,7 @@ void GenerateScrewDislocation(InArgs_t *inArgs)

		NormalizeVec(gDir);
		NormalizeVec(line);
		cross(gDir, line, normal);
		cross(line, gDir, normal);
		NormalizeVec(normal);
		FormatVector(normal, "normal");

		@@ -81,7 +81,11 @@ void GenerateScrewDislocation(InArgs_t *inArgs)

		if(inArgs->help)return;

		if(strstr(inArgs->inpFiles[0].c_str(), ".lmp") != NULL){
		ReadLMPFile(inArgs->inpFiles[0], dum);
		}else{
		ReadDumpFile(inArgs->inpFiles[0], dum);
		}

		for(i=0; i<3; i++){
		boundMin[i] = dum.box[i][0];
		@@ -90,6 +94,7 @@ void GenerateScrewDislocation(InArgs_t *inArgs)

		real8 point[3], dvec[3], dis, vec[3], a, b, shift;

		if(nPairs > 0){
		while(nPairs > 0){

		#pragma omp parallel for private(point, dvec, dis, vec, a, b, shift) shared(dum)
		@@ -121,7 +126,7 @@ void GenerateScrewDislocation(InArgs_t *inArgs)
		dum.atom[i].y += (shift*line[1]);
		dum.atom[i].z += (shift*line[2]);

		// FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
		FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
		}

		printf("pair %d: burgMag %f, position {%f,%f,%f}, rsize %f\n", nPairs, burgMag, pos[0], pos[1], pos[2], rsize);
		@@ -133,15 +138,165 @@ void GenerateScrewDislocation(InArgs_t *inArgs)

		rsize = delta;
		nPairs--;

		}
		}else{
		for(i=0; i<dum.atom.size(); i++){
		FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
		}
		}

		WriteDumpFile(inArgs->outFiles[0], dum);
		// WriteDumpFile(inArgs->outFiles[0], dum);
		MGToLMPDataFile(inArgs->outFiles[0], dum);
		return;
		}


		#define DisDisPlacement(r) (0.5 + sqrt(M_PI)(r)(15.0+10.0M_PI(r)(r) + 2.0M_PIM_PI(r)(r)(r)*(r)) \
		/ (4.0pow((2.0 + M_PI(r)*(r)), 2.5)) )

		void GenerateExtendedDislocation(InArgs_t *inArgs)
		{
		int index, nDiss, i;
		string posName("pos"), lineName("line"), separationName("separation"), gDirName("gdir");
		string burgName("burg");
		real8 boundMax[3], boundMin[3], pos[3];
		real8 line[3] = {1,0,0}, separation = 10.0, gDir[3] = {0,1,0}, normal[3];
		vector<vector<real8> > burgList, posList;

		if((index = GetValID(inArgs->priVars, posName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() == 3){
		pos[0] = atof(inArgs->priVars[index].vals[0].c_str());
		pos[1] = atof(inArgs->priVars[index].vals[1].c_str());
		pos[2] = atof(inArgs->priVars[index].vals[2].c_str());
		}
		}
		printf("The drag dislocation position (pos): {%f, %f, %f}\n", pos[0], pos[1], pos[2]);

		if((index = GetValID(inArgs->priVars, separationName)) < inArgs->priVars.size()){
		separation = atof(inArgs->priVars[index].vals[0].c_str());
		}
		printf("The separation of dislocations (separation): %f\n", separation);

		if((index = GetValID(inArgs->priVars, lineName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() == 3){
		line[0] = atof(inArgs->priVars[index].vals[0].c_str());
		line[1] = atof(inArgs->priVars[index].vals[1].c_str());
		line[2] = atof(inArgs->priVars[index].vals[2].c_str());
		}
		NormalizeVec(line);
		}
		printf("The line direction of dislocation (line): {%f, %f, %f}\n", line[0], line[1], line[2]);

		if((index = GetValID(inArgs->priVars, gDirName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() == 3){
		gDir[0] = atof(inArgs->priVars[index].vals[0].c_str());
		gDir[1] = atof(inArgs->priVars[index].vals[1].c_str());
		gDir[2] = atof(inArgs->priVars[index].vals[2].c_str());
		}
		NormalizeVec(gDir);
		}
		printf("The glide direction of dislocation (gDir): {%f, %f, %f}\n", gDir[0], gDir[1], gDir[2]);

		if((index = GetValID(inArgs->priVars, burgName)) < inArgs->priVars.size()){

		nDiss = (int)(inArgs->priVars[index].vals.size());
		nDiss /= 3;
		burgList.resize(nDiss);
		posList.resize(nDiss);

		printf("%d dislocations: \n", nDiss);
		for(i=0; i<burgList.size(); i++){
		burgList[i].resize(3);
		posList[i].resize(3);
		burgList[i][0] = atof(inArgs->priVars[index].vals[3*i].c_str());
		burgList[i][1] = atof(inArgs->priVars[index].vals[3*i+1].c_str());
		burgList[i][2] = atof(inArgs->priVars[index].vals[3*i+2].c_str());

		posList[i][0] = pos[0] + 2gDir[0]separation*((real8)i);
		posList[i][1] = pos[1] + 2gDir[1]separation*((real8)i);
		posList[i][2] = pos[2] + 2gDir[2]separation*((real8)i);

		printf("%d: position {%f,%f,%f}, burg {%f,%f,%f}\n", i, posList[i][0], posList[i][1], posList[i][2],
		burgList[i][0], burgList[i][1], burgList[i][2]);
		}
		}else{
		Fatal("use -dburg to define dislocation lists");
		}

		if(inArgs->help)return;
		printf("Reading atoms file %s ...\n", inArgs->inpFiles[0].c_str());

		Dump_t dum;
		if(strstr(inArgs->inpFiles[0].c_str(), ".lmp") != NULL){
		ReadLMPFile(inArgs->inpFiles[0], dum);
		}else{
		ReadDumpFile(inArgs->inpFiles[0], dum);
		}

		for(i=0; i<3; i++){
		boundMin[i] = dum.box[i][0];
		boundMax[i] = dum.box[i][1];
		}
		printf("%d atoms have been read\n", (int)dum.atom.size());

		if(fabs(DotProduct(line, gDir)) > 1.0E-2){
		Fatal("glide direction should be perpendicular to line direction.");
		}

		cross(line, gDir, normal);
		NormalizeVec(normal);
		printf("The glide plane is {%f,%f,%f}\n", normal[0], normal[1], normal[2]);

		int j;
		real8 point[3], vec1[3], vec2[3], displace, vec[3], dis1, dis2;
		real8 a, b, angle;


		#pragma omp parallel for shared(dum) private (j,point,vec,vec1,vec2,dis1,dis2,displace,a,b, angle)
		for(i=0; i<dum.atom.size(); i++){

		point[0] = dum.atom[i].x;
		point[1] = dum.atom[i].y;
		point[2] = dum.atom[i].z;

		for(j=0; j<nDiss; j++){
		VECTOR_COPY(vec, posList[j]);

		PointLineIntersection(point, line, vec, vec1, &dis1);
		vec1[0] -= point[0];
		vec1[1] -= point[1];
		vec1[2] -= point[2];

		NormalizeVec(vec1);

		a = DotProduct(vec1, gDir);
		b = DotProduct(vec1, normal);

		angle = atan2(b, a)/M_PI;
		if(angle < 0)angle += 1;
		angle /= 2.0;

		displace = angle*DisDisPlacement(dis1);
		if(isnan(displace)){
		displace = 0;
		printf("0");
		}
		dum.atom[i].x += burgList[j][0]*displace;
		dum.atom[i].y += burgList[j][1]*displace;
		dum.atom[i].z += burgList[j][2]*displace;

		}
		// FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
		}

		MGToLMPDataFile(inArgs->outFiles[0], dum);

		if(inArgs->outFiles.size() == 2){
		WriteDumpFile(inArgs->outFiles[1], dum);
		}
		return;
		}

		void GenerateDislocation(InArgs_t *inArgs){
		int index;
		string disTypeName("type");
		@@ -158,8 +313,12 @@ void GenerateDislocation(InArgs_t *inArgs){
		GenerateScrewDislocation(inArgs);
		break;

		case 1:
		GenerateExtendedDislocation(inArgs);
		break;

		default:
		GenerateScrewDislocation(inArgs);
		GenerateExtendedDislocation(inArgs);
		break;
		}

src/HandleExtendedDislocation.cpp

+61 −18

Original line number	Diff line number	Diff line
		@@ -392,6 +392,8 @@ typedef struct {
		typedef struct {
		long int timestep;
		double x,y,z;
		double absy;
		double vy;
		double separation;
		vector<double> vars;
		}EDState_t;
		@@ -410,12 +412,16 @@ bool com(Atom_t p, Atom_t q){

		bool com2(Probe_t p, Probe_t q)
		{
		return(p.y<p.y);
		return(p.y<q.y);
		}
		bool com3(EDState_t p, EDState_t q)
		{
		return(p.timestep<q.timestep);
		}
		bool com4(Probe_t p, Probe_t q)
		{
		return(p.nbr.size() > q.nbr.size());
		}

		void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		{
		@@ -517,21 +523,21 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		}

		for(i=0; i<dum.atom.size(); ){
		if(dum.atom[i].vars[indexVar] != dval){
		if(dum.atom[i].x < dum.box[0][0] + 5 \|\|
		dum.atom[i].y < dum.box[1][0] + 5 \|\|
		dum.atom[i].z < dum.box[2][0] + 5 \|\|
		dum.atom[i].x > dum.box[0][1] - 5 \|\|
		dum.atom[i].y > dum.box[1][1] - 5 \|\|
		dum.atom[i].z > dum.box[2][1] - 5){
		dum.atom.erase(dum.atom.begin()+i);
		continue;
		}

		if(dum.atom[i].x < dum.box[0][0] + 2 \|\|
		dum.atom[i].y < dum.box[1][0] + 2 \|\|
		dum.atom[i].z < dum.box[2][0] + 2 \|\|
		dum.atom[i].x > dum.box[0][1] - 2 \|\|
		dum.atom[i].y > dum.box[1][1] - 2 \|\|
		dum.atom[i].z > dum.box[2][1] - 2){
		if(dum.atom[i].vars[indexVar] != dval){
		dum.atom.erase(dum.atom.begin()+i);
		continue;
		}


		i++;
		}

		@@ -569,7 +575,7 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		lastIndex = i;
		firstNbr = 0;
		}
		probe.nbr.push_back(&dum.atom[i]);
		probe.nbr.push_back(&(dum.atom[i]));
		}
		}
		if(probe.nbr.size() > nums[0] && probe.nbr.size()<nums[1]){
		@@ -603,12 +609,13 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		probes[i].z /= (double)probes[i].nbr.size();

		}
		sort(probes.begin(), probes.end(), com2);

		double separation = 0;

		if(probes.size() > 0){
		separation = fabs(probes[0].y-probes.back().y);
		if(effeSepRange[0] > 0 && effeSepRange[1] >0 ){
		sort(probes.begin(), probes.end(), com2);
		if(separation > effeSepRange[0] && separation < effeSepRange[1] ){
		states[file].x = (probes[0].x+probes.back().x)*0.5;
		states[file].y = (probes[0].y+probes.back().y)*0.5;
		@@ -626,13 +633,33 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		}
		}
		}
		}else{
		sort(probes.begin(), probes.end(), com4);

		states[file].x = probes[0].x;
		states[file].y = probes[0].y;
		states[file].z = probes[0].z;
		states[file].timestep = dum.timestep;
		states[file].separation = 0;
		printf("file %d, nbrsize %d, timestep %d %f \n",file,(int)probes[0].nbr.size(), states[file].timestep, states[file].y);
		if(logfile ==1){
		states[file].vars.resize(list.variables.size()-1);
		for(i=0; i<list.data[0].size(); i++){
		if(states[file].timestep == ((int)list.data[0][i]))break;
		}
		for(j=0; j<states[file].vars.size(); j++){
		states[file].vars[j] = list.data[j+1][i];
		}
		}

		}
		}

		}

		sort(states.begin(), states.end(), com3);

		out << "variables = timestep, separation, y";
		out << "variables = timestep, separation, y, absy, vy";
		if(logfile == 1){
		for(i=1; i<list.variables.size(); i++){
		out << ", " << list.variables[i];
		@@ -640,12 +667,28 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		}
		out << endl;

		real8 lastAbsy, vy, lasty, lastTimestep;
		ReadDumpFile(inArgs->inpFiles[0], dum);
		real8 ybox = dum.box[1][1];
		printf("effective states %d, boxy is %f\n", (int)states.size(), dum.box[1][1]);
		for(i=0; i<states.size(); i++){
		if(states[i].timestep > 1E9)break;

		lastAbsy = ((i==0) ? states[i].y : states[i-1].absy);
		lasty = ((i==0) ? states[i].y : states[i-1].y);
		lastTimestep = ((i==0) ? states[i].timestep : states[i-1].timestep);
		out << states[i].timestep << " ";
		out << states[i].separation << " ";
		out << states[i].y << " ";

		vy = states[i].y - lasty;
		vy -= rint(vy * 1.0/ybox) * ybox;
		states[i].absy = lastAbsy+vy;
		out << states[i].absy << " ";

		states[i].vy = ((i==0) ? 0 : (vy)/(states[i].timestep - lastTimestep));
		out << states[i].vy << " ";

		for(j=0; j< states[i].vars.size(); j++){
		out << setprecision(10) << states[i].vars[j] << " ";
		}

src/Math.cpp

+6 −9

Original line number	Diff line number	Diff line
		@@ -410,18 +410,20 @@ int PointLineIntersection(double p1, double d, double *p2,

		norm2 = sqrt(SQUARE(p1x-p2x) + SQUARE(p1y-p2y) + SQUARE(p1z-p2z));
		if(norm2 < EPS1){
		vec[0] = p1[0];
		vec[1] = p1[1];
		vec[2] = p1[2];
		vec[0] = p2[0];
		vec[1] = p2[1];
		vec[2] = p2[2];
		*dis = 0.0;
		printf("0");
		return(POINT_INTERSECTION);
		}

		if(fabs(dx(p1x-p2x) + dy(p1y-p2y) + dz(p1z-p2z)) < EPS1 norm *norm2){
		if(fabs(dx(p1x-p2x) + dy(p1y-p2y) + dz(p1z-p2z)) > normnorm2*(1-EPS1)){
		vec[0] = p1[0];
		vec[1] = p1[1];
		vec[2] = p1[2];
		*dis = 0.0;
		printf("0");
		return(POINT_INTERSECTION);
		}

		@@ -435,14 +437,9 @@ int PointLineIntersection(double p1, double d, double *p2,
		if(isnan(dis) \|\| isinf(dis)){
		printf(" PointLineIntersection: %f {%f,%f,%f}\n",
		*dis, vec[0], vec[1], vec[2]);

		FormatVector(p1, "p1");
		FormatVector(d, "d");
		FormatVector(p2, "p2");
		vec[0] = p1[0];
		vec[1] = p1[1];
		vec[2] = p1[2];
		Fatal("\n");
		}

		if(isnan(dis))dis = 0.0;