MDStressLab++
Mls.cpp
Go to the documentation of this file.
1 /*
2  * Mls.cpp
3  *
4  * Created on: Feb 11, 2020
5  * Author: Min
6  */
7 
8 #include <fstream>
9 #include <vector>
10 #include <math.h>
11 #include "Grid.h"
12 #include "BoxConfiguration.h"
13 #include "Configuration.h"
14 #include "SubConfiguration.h"
15 #include "Mls.h"
16 #include "neighbor_list.h"
17 #include "typedef.h"
18 
19 //Mls::Mls(const BoxConfiguration& body, const std::vector<Vector3d>& gridCoordinates, \
20  double radiusMls, const std::string name):radiusMls(radiusMls),name(name)
21 Mls::Mls(const BoxConfiguration& body, const Grid<Reference>* pgrid, \
22  double radiusMls, const std::string name):radiusMls(radiusMls),name(name)
23 {
24  MY_BANNER("Constructing Deformation Gradient using the Moving Least Square method!");
25  deformationGradient.reserve(pgrid->coordinates.size());
26  gridPushed.reserve(pgrid->coordinates.size());
27 
28  std::unique_ptr<const Configuration> pconfigMls;
29  if (body.pbc.any() == 1)
30  {
31  MY_HEADING("Generating padding atoms for periodic boundary conditions of Moving Least Squares.")
32  pconfigMls.reset(body.getConfiguration(radiusMls));
33  std::cout << "Thickness of padding = radiusMls = "
34  << radiusMls << std::endl;
35  std::cout << "Total number of atoms including padding atoms = " << pconfigMls->numberOfParticles << std::endl;
36  std::cout << std::endl;
37  }
38  else
39  {
40  pconfigMls.reset(body.getConfiguration(0.0));
41  //pconfigMls.reset(&body);
42  }
43  double box_xx, box_yy, box_zz;
44  box_xx = body.box(0);
45  box_yy = body.box(4);
46  box_zz = body.box(8);
47  int pbc_x, pbc_y, pbc_z;
48  pbc_x = body.pbc(0);
49  pbc_y = body.pbc(1);
50  pbc_z = body.pbc(2);
51  //MatrixXd displacements(body.coordinates.at(Reference).rows(),body.coordinates.at(Reference).cols());
52  MatrixXd displacements(pconfigMls->coordinates.at(Reference).rows(),pconfigMls->coordinates.at(Reference).cols());
53  displacements.setZero();
54 /*
55  // need to respect boundary conditions
56  for (int i = 0; i != body.coordinates.at(Reference).rows(); i++)
57  {
58 
59  if (fabs(body.coordinates.at(Current)(i,0) - body.coordinates.at(Reference)(i,0)) > 0.5 * body.box(0) && body.pbc(0))
60  {
61  displacements(i,0) = body.coordinates.at(Current)(i,0) - body.coordinates.at(Reference)(i,0) - body.box(0);
62  }
63  else
64  {
65  displacements(i,0) = body.coordinates.at(Current)(i,0) - body.coordinates.at(Reference)(i,0);
66  }
67 
68  if (fabs(body.coordinates.at(Current)(i,1) - body.coordinates.at(Reference)(i,1)) > 0.5 * body.box(4) && body.pbc(1))
69  {
70  displacements(i,1) = body.coordinates.at(Current)(i,1) - body.coordinates.at(Reference)(i,1) - body.box(4);
71  }
72  else
73  {
74  displacements(i,1) = body.coordinates.at(Current)(i,1) - body.coordinates.at(Reference)(i,1);
75  }
76 
77  if (fabs(body.coordinates.at(Current)(i,2) - body.coordinates.at(Reference)(i,2)) > 0.5 * body.box(8) && body.pbc(2))
78  {
79  displacements(i,2) = body.coordinates.at(Current)(i,2) - body.coordinates.at(Reference)(i,2) - body.box(8);
80  }
81  else
82  {
83  displacements(i,2) = body.coordinates.at(Current)(i,2) - body.coordinates.at(Reference)(i,2);
84  }
85  }
86 
87  //std::cout << "Box size " << body.box(0) << " " << body.box(4) << " " << body.box(8) << " " << std::endl;
88 */
89 /*
90  // need to respect boundary conditions
91  for (int i = 0; i != pconfigMls->coordinates.at(Reference).rows(); i++)
92  {
93 
94  if (fabs(pconfigMls->coordinates.at(Current)(i,0) - pconfigMls->coordinates.at(Reference)(i,0)) > 0.5 * box_xx && pbc_x)
95  {
96  displacements(i,0) = pconfigMls->coordinates.at(Current)(i,0) - pconfigMls->coordinates.at(Reference)(i,0) - box_xx;
97  }
98  else
99  {
100  displacements(i,0) = pconfigMls->coordinates.at(Current)(i,0) - pconfigMls->coordinates.at(Reference)(i,0);
101  }
102 
103  if (fabs(pconfigMls->coordinates.at(Current)(i,1) - pconfigMls->coordinates.at(Reference)(i,1)) > 0.5 * box_yy && pbc_y)
104  {
105  displacements(i,1) = pconfigMls->coordinates.at(Current)(i,1) - pconfigMls->coordinates.at(Reference)(i,1) - box_yy;
106  }
107  else
108  {
109  displacements(i,1) = pconfigMls->coordinates.at(Current)(i,1) - pconfigMls->coordinates.at(Reference)(i,1);
110  }
111 
112  if (fabs(pconfigMls->coordinates.at(Current)(i,2) - pconfigMls->coordinates.at(Reference)(i,2)) > 0.5 * box_zz && pbc_z)
113  {
114  displacements(i,2) = pconfigMls->coordinates.at(Current)(i,2) - pconfigMls->coordinates.at(Reference)(i,2) - box_zz;
115  }
116  else
117  {
118  displacements(i,2) = pconfigMls->coordinates.at(Current)(i,2) - pconfigMls->coordinates.at(Reference)(i,2);
119  }
120  }
121 */
122 
123  double gridRadiusMls;
124  double r2;
125  Vector3d r;
126  Vector3d rSanityI, rSanityJ, rSanityK, rSanityQ;
127  int sanityI, sanityJ, sanityK, sanityQ;
128  bool sanityCheck;
129  int cycle_count;
130 
131  Matrix3d tensorF;
132  Vector3d gptIPushedF;
133  Matrix4d sanity, A, inverseA, dAdx, dAdy, dAdz;
134  Matrix4d inverseAXdAdx, inverseAXdAdy, inverseAXdAdz;
135  Matrix4d inverseAXdAdxXInverseA, inverseAXdAdyXInverseA, inverseAXdAdzXInverseA;
136  MatrixXd MlsDisp(4,3), dMlsDispdx(4,3), dMlsDispdy(4,3), dMlsDispdz(4,3);
137  MatrixXd fintmdx, fintmdy,fintmdz, sintmdx, sintmdy, sintmdz;
138  MatrixXd B, inverseAXB, dBdx, dBdy, dBdz;
139  MatrixXd exactDisplacements;
140  MatrixXd P, dWdx;
141  VectorXd W;
142 
143  std::vector<int> gridMlsAtomList;
144  int ngridMLSAtomList;
145 
146  // ------------------------------------------------------------------
147  // TODO: Trying to build grid neighbor list for MLS
148  // ------------------------------------------------------------------
149  Stencil stencil(*pconfigMls);
150  //stencil.expandStencil(pgrid,10.0 * radiusMls,10.0 * radiusMls);
151  stencil.expandStencil(pgrid,radiusMls,0.0);
152  SubConfiguration subconfig{stencil};
153  //std::cout << "debug subconfig number of particles: " << subconfig.numberOfParticles << std::endl;
154 
155 /*
156  for (auto& [key, value]: subconfig.globalLocalMap)
157  {
158  std::cout << "global paricle id =" << key+1 << ", local particle id = " << value+1 << std::endl;
159  }
160 */
161 
162  std::vector<std::set<int>> neighborListsOfGridsMls=pgrid->getGridNeighborLists(subconfig,radiusMls);
163 
164  // need to respect boundary conditions
165  for (int i = 0; i != subconfig.coordinates.at(Reference).rows(); i++)
166  {
167 
168  if (fabs(subconfig.coordinates.at(Current)(i,0) - subconfig.coordinates.at(Reference)(i,0)) > 0.5 * box_xx && pbc_x)
169  {
170  displacements(i,0) = subconfig.coordinates.at(Current)(i,0) - subconfig.coordinates.at(Reference)(i,0) - box_xx;
171  }
172  else
173  {
174  displacements(i,0) = subconfig.coordinates.at(Current)(i,0) - subconfig.coordinates.at(Reference)(i,0);
175  }
176 
177  if (fabs(subconfig.coordinates.at(Current)(i,1) - subconfig.coordinates.at(Reference)(i,1)) > 0.5 * box_yy && pbc_y)
178  {
179  displacements(i,1) = subconfig.coordinates.at(Current)(i,1) - subconfig.coordinates.at(Reference)(i,1) - box_yy;
180  }
181  else
182  {
183  displacements(i,1) = subconfig.coordinates.at(Current)(i,1) - subconfig.coordinates.at(Reference)(i,1);
184  }
185 
186  if (fabs(subconfig.coordinates.at(Current)(i,2) - subconfig.coordinates.at(Reference)(i,2)) > 0.5 * box_zz && pbc_z)
187  {
188  displacements(i,2) = subconfig.coordinates.at(Current)(i,2) - subconfig.coordinates.at(Reference)(i,2) - box_zz;
189  }
190  else
191  {
192  displacements(i,2) = subconfig.coordinates.at(Current)(i,2) - subconfig.coordinates.at(Reference)(i,2);
193  }
194  }
195 
196  //std::cout << "debug: number of atoms in the neighborlist: " << neighborListsOfGridsMls[0].size() << std::endl;
197 
198  //for (std::vector<Vector3d>::size_type iGrid = 0; iGrid != gridCoordinates.size(); iGrid++)
199  for (int iGrid = 0; iGrid != pgrid->coordinates.size(); iGrid++)
200  //for (std::vector<Vector3d>::size_type iGrid = 0; iGrid != 1; iGrid++) // debug only calculate one point
201  {
202  gridRadiusMls = radiusMls;
203  A = Matrix4d::Zero();
204  cycle_count = 1;
205 
206  do
207  {
208  if (cycle_count != 1)
209  {
210  //MY_WARNING("Something is wrong and need to double the MLS_Radius for grid point " + std::to_string(iGrid + 1) \
211  + ". Starting Cycle " + std::to_string(cycle_count) + ".")
212  MY_WARNING("Something is wrong and causing the A matrix singular, \
213  Probably because there are not enough atoms in the MLS radius. Setting the deformation gradient to be identity and \
214  pushed grid point to the grid point itself.")
215  tensorF.setIdentity();
216  /*
217  tensorF(0,0) = 1.0;
218  tensorF(0,1) = 0.0;
219  tensorF(0,2) = 0.0;
220  tensorF(1,0) = 0.0;
221  tensorF(1,1) = 1.0;
222  tensorF(1,2) = 0.0;
223  tensorF(2,0) = 0.0;
224  tensorF(2,1) = 0.0;
225  tensorF(2,2) = 1.0 ;
226  */
227  gptIPushedF = pgrid->coordinates[iGrid];
228 
229  deformationGradient.push_back(tensorF.transpose());
230  gridPushed.push_back(gptIPushedF);
231  goto GridEnd;
232  }
233 
234  //if (ngridMLSAtomList == body.coordinates.at(Reference).rows())
235  if (ngridMLSAtomList == pconfigMls->coordinates.at(Reference).rows())
236  {
237  MY_ERROR("Model degenerate to 2D. Cannot use MLS3D to generate deformation gradient!")
238  }
239 
240  gridMlsAtomList.clear();
241  ngridMLSAtomList = 0;
242  const std::set<int>& neighborListOfGridsMls = neighborListsOfGridsMls[iGrid];
243  //for (int jRefAtoms = 0; jRefAtoms != body.coordinates.at(Reference).rows(); jRefAtoms++)
244  //for (int jRefAtoms = 0; jRefAtoms != pconfigMls->coordinates.at(Reference).rows(); jRefAtoms++)
245  //for (auto jRefAtoms = neighborListsOfGridsMls[iGrid].cbegin(); jRefAtoms != neighborListsOfGridsMls[iGrid].cend(); jRefAtoms++)
246  for (const auto& jRefAtoms : neighborListOfGridsMls)
247  {
248  //r(0) = body.coordinates.at(Reference)(jRefAtoms,0) - gridCoordinates[iGrid](0);
249  //r(1) = body.coordinates.at(Reference)(jRefAtoms,1) - gridCoordinates[iGrid](1);
250  //r(2) = body.coordinates.at(Reference)(jRefAtoms,2) - gridCoordinates[iGrid](2);
251  r(0) = subconfig.coordinates.at(Reference)(jRefAtoms,0) - pgrid->coordinates[iGrid](0);
252  r(1) = subconfig.coordinates.at(Reference)(jRefAtoms,1) - pgrid->coordinates[iGrid](1);
253  r(2) = subconfig.coordinates.at(Reference)(jRefAtoms,2) - pgrid->coordinates[iGrid](2);
254  //r(0) = pconfigMls->coordinates.at(Reference)(jRefAtoms,0) - pgrid->coordinates[iGrid](0);
255  //r(1) = pconfigMls->coordinates.at(Reference)(jRefAtoms,1) - pgrid->coordinates[iGrid](1);
256  //r(2) = pconfigMls->coordinates.at(Reference)(jRefAtoms,2) - pgrid->coordinates[iGrid](2);
257  //std::cout <<"debug push_back1: " << jRefAtoms+1 << std::endl;
258 
259  if (r.norm() <= gridRadiusMls)
260  {
261 
262  //std::cout << "debug: positions: " << r << std::endl;
263  //std::cout << r.norm() << std::endl;
264  //gridMlsAtomList.push_back(jRefAtoms);
265  //std::cout <<"debug push_back: " << jRefAtoms << std::endl;
266  gridMlsAtomList.push_back(jRefAtoms);
267  //std::cout <<"debug push_back: " << jRefAtoms+1 << std::endl;
268  ngridMLSAtomList++;
269  }
270  }
271  //exit(0);
272  // sanity check start!
273  sanityCheck = false;
274 
275  if (ngridMLSAtomList <= 3)
276  {
277  sanityCheck = false;
278  A = Matrix4d::Zero();
279  MY_WARNING("Not enough atoms in the MLS radius " + std::to_string(gridRadiusMls) + ". Doubling it.");
280  gridRadiusMls = 2.0 * gridRadiusMls;
281  goto EndLoop;
282  }
283  else
284  {
285  for (int i = 0; i <= ngridMLSAtomList - 4; i++)
286  {
287  sanityI = gridMlsAtomList[i];
288  //rSanityI(0) = body.coordinates.at(Reference)(sanityI,0) - gridCoordinates[iGrid](0);
289  //rSanityI(1) = body.coordinates.at(Reference)(sanityI,1) - gridCoordinates[iGrid](1);
290  //rSanityI(2) = body.coordinates.at(Reference)(sanityI,2) - gridCoordinates[iGrid](2);
291  rSanityI(0) = subconfig.coordinates.at(Reference)(sanityI,0) - pgrid->coordinates[iGrid](0);
292  rSanityI(1) = subconfig.coordinates.at(Reference)(sanityI,1) - pgrid->coordinates[iGrid](1);
293  rSanityI(2) = subconfig.coordinates.at(Reference)(sanityI,2) - pgrid->coordinates[iGrid](2);
294  for(int j = i + 1; j <= ngridMLSAtomList - 3; j++)
295  {
296  sanityJ = gridMlsAtomList[j];
297  //rSanityJ(0) = body.coordinates.at(Reference)(sanityJ,0) - gridCoordinates[iGrid](0);
298  //rSanityJ(1) = body.coordinates.at(Reference)(sanityJ,1) - gridCoordinates[iGrid](1);
299  //rSanityJ(2) = body.coordinates.at(Reference)(sanityJ,2) - gridCoordinates[iGrid](2);
300  rSanityJ(0) = subconfig.coordinates.at(Reference)(sanityJ,0) - pgrid->coordinates[iGrid](0);
301  rSanityJ(1) = subconfig.coordinates.at(Reference)(sanityJ,1) - pgrid->coordinates[iGrid](1);
302  rSanityJ(2) = subconfig.coordinates.at(Reference)(sanityJ,2) - pgrid->coordinates[iGrid](2);
303  for (int k = j + 1; k <= ngridMLSAtomList - 2; k++)
304  {
305  sanityK = gridMlsAtomList[k];
306  //rSanityK(0) = body.coordinates.at(Reference)(sanityK,0) - gridCoordinates[iGrid](0);
307  //rSanityK(1) = body.coordinates.at(Reference)(sanityK,1) - gridCoordinates[iGrid](1);
308  //rSanityK(2) = body.coordinates.at(Reference)(sanityK,2) - gridCoordinates[iGrid](2);
309  rSanityK(0) = subconfig.coordinates.at(Reference)(sanityK,0) - pgrid->coordinates[iGrid](0);
310  rSanityK(1) = subconfig.coordinates.at(Reference)(sanityK,1) - pgrid->coordinates[iGrid](1);
311  rSanityK(2) = subconfig.coordinates.at(Reference)(sanityK,2) - pgrid->coordinates[iGrid](2);
312  for (int q = k + 1; q <= ngridMLSAtomList - 1; q++)
313  {
314  sanityQ = gridMlsAtomList[q];
315  //rSanityQ(0) = body.coordinates.at(Reference)(sanityQ,0) - gridCoordinates[iGrid](0);
316  //rSanityQ(1) = body.coordinates.at(Reference)(sanityQ,1) - gridCoordinates[iGrid](1);
317  //rSanityQ(2) = body.coordinates.at(Reference)(sanityQ,2) - gridCoordinates[iGrid](2);
318  rSanityQ(0) = subconfig.coordinates.at(Reference)(sanityQ,0) - pgrid->coordinates[iGrid](0);
319  rSanityQ(1) = subconfig.coordinates.at(Reference)(sanityQ,1) - pgrid->coordinates[iGrid](1);
320  rSanityQ(2) = subconfig.coordinates.at(Reference)(sanityQ,2) - pgrid->coordinates[iGrid](2);
321 
322 
323  sanity = Matrix4d::Constant(1.0);
324  sanity(0,0) = rSanityI(0);
325  sanity(0,1) = rSanityI(1);
326  sanity(0,2) = rSanityI(2);
327  sanity(1,0) = rSanityJ(0);
328  sanity(1,1) = rSanityJ(1);
329  sanity(1,2) = rSanityJ(2);
330  sanity(2,0) = rSanityK(0);
331  sanity(2,1) = rSanityK(1);
332  sanity(2,2) = rSanityK(2);
333  sanity(3,0) = rSanityQ(0);
334  sanity(3,1) = rSanityQ(1);
335  sanity(3,2) = rSanityQ(2);
336 
337  if (sanity.determinant() > epsilon)
338  {
339  sanityCheck = true;
340  //inverseA = Matrix4d::Identity();
341  //std::cout << "Passed Sanity Check!" << std::endl;
342  goto SanityDone;
343  }
344  }
345  }
346  }
347  }
348  if (!sanityCheck)
349  {
350  A = Matrix4d::Zero();
351  MY_WARNING("All atoms lying on a surface for the MLS radius " + std::to_string(gridRadiusMls) + ". Doubling it.");
352  gridRadiusMls = 2.0 * gridRadiusMls;
353  goto EndLoop;
354  }
355  }
356 
357  SanityDone:;
358 
359  sort(gridMlsAtomList.begin(), gridMlsAtomList.end());
360 
361  //for (auto it = gridMlsAtomList.begin(); it != gridMlsAtomList.end(); ++it)
362  //{
363  // std::cout << ' ' << *it + 1 << ' ';
364  //}
365  //std::cout << std::endl;
366  //std::cout << "ngridMLSAtomList: " << ngridMLSAtomList << std::endl;
367 
368  P.resize(ngridMLSAtomList,4);
369  W.resize(ngridMLSAtomList);
370  dWdx.resize(ngridMLSAtomList,3);
371  B.resize(4,ngridMLSAtomList);
372  dBdx.resize(4,ngridMLSAtomList);
373  dBdy.resize(4,ngridMLSAtomList);
374  dBdz.resize(4,ngridMLSAtomList);
375  inverseAXB.resize(4,ngridMLSAtomList);
376  exactDisplacements.resize(ngridMLSAtomList,3);
377  fintmdx.resize(4,ngridMLSAtomList);
378  fintmdy.resize(4,ngridMLSAtomList);
379  fintmdz.resize(4,ngridMLSAtomList);
380  sintmdx.resize(4,ngridMLSAtomList);
381  sintmdy.resize(4,ngridMLSAtomList);
382  sintmdz.resize(4,ngridMLSAtomList);
383 
384  tensorF = Matrix3d::Zero();
385  gptIPushedF = Vector3d::Zero();
386  P = MatrixXd::Zero(ngridMLSAtomList,4);
387  W = VectorXd::Zero(ngridMLSAtomList);
388  dWdx = MatrixXd::Zero(ngridMLSAtomList,3);
389  A = Matrix4d::Zero();
390  inverseA = Matrix4d::Zero();
391  dAdx = Matrix4d::Zero();
392  dAdy = Matrix4d::Zero();
393  dAdz = Matrix4d::Zero();
394  B = MatrixXd::Zero(4,ngridMLSAtomList);
395  dBdx = MatrixXd::Zero(4,ngridMLSAtomList);
396  dBdy = MatrixXd::Zero(4,ngridMLSAtomList);
397  dBdz = MatrixXd::Zero(4,ngridMLSAtomList);
398  inverseAXB = MatrixXd::Zero(4,ngridMLSAtomList);
399  MlsDisp = MatrixXd::Zero(4,3);
400  dMlsDispdx = MatrixXd::Zero(4,3);
401  dMlsDispdy = MatrixXd::Zero(4,3);
402  dMlsDispdz = MatrixXd::Zero(4,3);
403  exactDisplacements = MatrixXd::Zero(ngridMLSAtomList,3);
404  inverseAXdAdx = Matrix4d::Zero();
405  inverseAXdAdy = Matrix4d::Zero();
406  inverseAXdAdz = Matrix4d::Zero();
407  inverseAXdAdxXInverseA = Matrix4d::Zero();
408  inverseAXdAdyXInverseA = Matrix4d::Zero();
409  inverseAXdAdzXInverseA = Matrix4d::Zero();
410  fintmdx = MatrixXd::Zero(4,ngridMLSAtomList);
411  fintmdy = MatrixXd::Zero(4,ngridMLSAtomList);
412  fintmdz = MatrixXd::Zero(4,ngridMLSAtomList);
413  sintmdx = MatrixXd::Zero(4,ngridMLSAtomList);
414  sintmdy = MatrixXd::Zero(4,ngridMLSAtomList);
415  sintmdz = MatrixXd::Zero(4,ngridMLSAtomList);
416 
417  for (int i = 0; i < ngridMLSAtomList; i++)
418  {
419  //r(0) = body.coordinates.at(Reference)(gridMlsAtomList[i],0) - gridCoordinates[iGrid](0);
420  //r(1) = body.coordinates.at(Reference)(gridMlsAtomList[i],1) - gridCoordinates[iGrid](1);
421  //r(2) = body.coordinates.at(Reference)(gridMlsAtomList[i],2) - gridCoordinates[iGrid](2);
422  r(0) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],0) - pgrid->coordinates[iGrid](0);
423  r(1) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],1) - pgrid->coordinates[iGrid](1);
424  r(2) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],2) - pgrid->coordinates[iGrid](2);
425  // Form P
426  P(i,0) = 1.0;
427  //P(i,1) = body.coordinates.at(Reference)(gridMlsAtomList[i],0) - gridCoordinates[iGrid](0);
428  //P(i,2) = body.coordinates.at(Reference)(gridMlsAtomList[i],1) - gridCoordinates[iGrid](1);
429  //P(i,3) = body.coordinates.at(Reference)(gridMlsAtomList[i],2) - gridCoordinates[iGrid](2);
430  P(i,1) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],0) - pgrid->coordinates[iGrid](0);
431  P(i,2) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],1) - pgrid->coordinates[iGrid](1);
432  P(i,3) = subconfig.coordinates.at(Reference)(gridMlsAtomList[i],2) - pgrid->coordinates[iGrid](2);
433  // Form W and dWdx
434  r2 = r.norm() / gridRadiusMls;
435  if (r2 <= 0.5)
436  {
437  W(i) = 2.0 / 3.0 - 4.0 * r2 * r2 + 4.0 * r2 * r2 * r2;
438  dWdx(i,0) = -8.0 * r(0) / gridRadiusMls / gridRadiusMls \
439  + 12.0 * r2 * r(0) / gridRadiusMls / gridRadiusMls;
440  dWdx(i,1) = -8.0 * r(1) / gridRadiusMls / gridRadiusMls \
441  + 12.0 * r2 * r(1) / gridRadiusMls / gridRadiusMls;
442  dWdx(i,2) = -8.0 * r(2) / gridRadiusMls / gridRadiusMls \
443  + 12.0 * r2 * r(2) / gridRadiusMls / gridRadiusMls;
444  }
445  else if (r2 >= 0.5 && r2 <= 1.0)
446  {
447  W(i) = 4.0 / 3.0 - 4.0 * r2 + 4.0 * r2 * r2 - 4.0 / 3.0 * r2 * r2 * r2;
448  dWdx(i,0) = -4.0 * r(0) / gridRadiusMls / gridRadiusMls / r2 \
449  + 8.0 * r(0) / gridRadiusMls / gridRadiusMls \
450  - 4.0 * r2 * r(0) / gridRadiusMls / gridRadiusMls;
451  dWdx(i,1) = -4.0 * r(1) / gridRadiusMls / gridRadiusMls / r2 \
452  + 8.0 * r(1) / gridRadiusMls / gridRadiusMls \
453  - 4.0 * r2 * r(1) / gridRadiusMls / gridRadiusMls;
454  dWdx(i,2) = -4.0 * r(2) / gridRadiusMls / gridRadiusMls / r2 \
455  + 8.0 * r(2) / gridRadiusMls / gridRadiusMls \
456  - 4.0 * r2 * r(2) / gridRadiusMls / gridRadiusMls;
457  }
458  else
459  {
460  W(i) = 0.0;
461  dWdx(i,0) = 0.0;
462  dWdx(i,1) = 0.0;
463  dWdx(i,2) = 0.0;
464  }
465  exactDisplacements(i,0) = displacements(gridMlsAtomList[i],0);
466  exactDisplacements(i,1) = displacements(gridMlsAtomList[i],1);
467  exactDisplacements(i,2) = displacements(gridMlsAtomList[i],2);
468  }
469 
470  for (int i = 0; i < 4; i++)
471  {
472  for (int j = 0; j < ngridMLSAtomList; j++)
473  {
474  B(i,j) = P(j,i) * W(j);
475  }
476  }
477 
478  for (int i = 0; i < 4; i++)
479  {
480  for (int j = 0; j < 4; j++)
481  {
482  for (int k = 0; k < ngridMLSAtomList; k++)
483  {
484  A(i,j) = A(i,j) + B(i,k) * P(k,j);
485  }
486  }
487  }
488 
489  inverseA = A.inverse();
490  gridRadiusMls = 2.0 * gridRadiusMls;
491 
492 /*
493  double detmat44;
494  detmat44 = A(1,1)*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))) \
495  - A(1,2)*(A(2,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))) \
496  + A(1,3)*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) \
497  - A(1,4)*(A(2,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(2,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(2,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)));
498 
499  inverseA(1,1) = \
500  (A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))) / detmat44;
501  inverseA(2,1) = \
502  (A(2,1)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(2,3)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(2,4)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))) / detmat44;
503  inverseA(3,1) = \
504  (A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) / detmat44;
505  inverseA(4,1) = \
506  (A(2,1)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))+A(2,2)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))+A(2,3)*(A(3,2)*A(4,1)-A(3,1)*A(4,2))) / detmat44;
507  inverseA(1,2) = \
508  (A(1,2)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(1,3)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(1,4)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))) / detmat44;
509  inverseA(2,2) = \
510  (A(1,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(1,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(1,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))) / detmat44;
511  inverseA(3,2) = \
512  (A(1,1)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(1,2)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(1,4)*(A(3,2)*A(4,1)-A(3,1)*A(4,2))) / detmat44;
513  inverseA(4,2) = \
514  (A(1,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(1,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(1,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) / detmat44;
515  inverseA(1,3) = \
516  (A(1,2)*(A(2,3)*A(4,4)-A(2,4)*A(4,3))+A(1,3)*(A(2,4)*A(4,2)-A(2,2)*A(4,4))+A(1,4)*(A(2,2)*A(4,3)-A(2,3)*A(4,2))) / detmat44;
517  inverseA(2,3) = \
518  (A(1,1)*(A(2,4)*A(4,3)-A(2,3)*A(4,4))+A(1,3)*(A(2,1)*A(4,4)-A(2,4)*A(4,1))+A(1,4)*(A(2,3)*A(4,1)-A(2,1)*A(4,3))) / detmat44;
519  inverseA(3,3) = \
520  (A(1,1)*(A(2,2)*A(4,4)-A(2,4)*A(4,2))+A(1,2)*(A(2,4)*A(4,1)-A(2,1)*A(4,4))+A(1,4)*(A(2,1)*A(4,2)-A(2,2)*A(4,1))) / detmat44;
521  inverseA(4,3) = \
522  (A(1,1)*(A(2,3)*A(4,2)-A(2,2)*A(4,3))+A(1,2)*(A(2,1)*A(4,3)-A(2,3)*A(4,1))+A(1,3)*(A(2,2)*A(4,1)-A(2,1)*A(4,2))) / detmat44;
523  inverseA(1,4) = \
524  (A(1,2)*(A(2,4)*A(3,3)-A(2,3)*A(3,4))+A(1,3)*(A(2,2)*A(3,4)-A(2,4)*A(3,2))+A(1,4)*(A(2,3)*A(3,2)-A(2,2)*A(3,3))) / detmat44;
525  inverseA(2,4) = \
526  (A(1,1)*(A(2,3)*A(3,4)-A(2,4)*A(3,3))+A(1,3)*(A(2,4)*A(3,1)-A(2,1)*A(3,4))+A(1,4)*(A(2,1)*A(3,3)-A(2,3)*A(3,1))) / detmat44;
527  inverseA(3,4) = \
528  (A(1,1)*(A(2,4)*A(3,2)-A(2,2)*A(3,4))+A(1,2)*(A(2,1)*A(3,4)-A(2,4)*A(3,1))+A(1,4)*(A(2,2)*A(3,1)-A(2,1)*A(3,2))) / detmat44;
529  inverseA(4,4) = \
530  (A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))) / detmat44;
531 
532 */
533 /*
534  //std::cout << std::endl << "ngridMLSAtomList: " << ngridMLSAtomList << std::endl;
535 
536  std::cout << "P" << std::endl;
537  std::cout << (P.transpose()).format(Eigen::FullPrecision) << std::endl;
538  std::cout << "W" << std::endl;
539  std::cout << W.format(Eigen::FullPrecision) << std::endl;
540 
541  std::cout << "B" << std::endl;
542  std::cout << B.transpose().format(Eigen::FullPrecision) << std::endl;
543 
544  std::cout << "exactDisplacements: " << std::endl;
545  std::cout << exactDisplacements.format(Eigen::FullPrecision) << std::endl;
546 
547  std::cout << "A: " << std::endl;
548  std::cout << (A.transpose()).format(Eigen::FullPrecision) << std::endl;
549 
550  std::cout << "A.determinant: " << A.determinant() << std::endl;
551 
552  std::cout << "inverseA: " << std::endl;
553  std::cout << (inverseA.transpose()).format(Eigen::FullPrecision) << std::endl;
554 
555  std::cout << "dWdx: " << std::endl;
556  std::cout << (dWdx.transpose()).format(Eigen::FullPrecision) << std::endl;
557 */
558 
559  EndLoop:;
560  cycle_count++;
561  } while (A.determinant() < epsilon);
562 
563  gridRadiusMls = gridRadiusMls / 2.0;
564 
565  // inv_A * B * exact_disp
566  for (int i = 0; i < 4; i++)
567  {
568  for (int j = 0; j < ngridMLSAtomList; j++)
569  {
570  for (int k = 0; k < 4; k++)
571  {
572  inverseAXB(i,j) = inverseAXB(i,j) + inverseA(i,k) * B(k,j);
573  }
574  }
575  }
576 
577  for (int i = 0; i < 4; i++)
578  {
579  for (int j = 0; j < 3; j++)
580  {
581  for (int k = 0; k < ngridMLSAtomList; k++)
582  {
583  MlsDisp(i,j) = MlsDisp(i,j) + inverseAXB(i,k) * exactDisplacements(k,j);
584  }
585  }
586  }
587 
588  //std::cout << "MlsDisp: " << std::endl;
589  //std::cout << (MlsDisp.transpose()).format(Eigen::FullPrecision) << std::endl;
590 
591  for (int i = 0; i < 4; i++)
592  {
593  for (int j = 0; j < ngridMLSAtomList; j++)
594  {
595  dBdx(i,j) = P(j,i) * dWdx(j,0);
596  dBdy(i,j) = P(j,i) * dWdx(j,1);
597  dBdz(i,j) = P(j,i) * dWdx(j,2);
598  }
599  }
600 
601  for (int i = 0; i < 4; i++)
602  {
603  for (int j = 0; j < 4; j++)
604  {
605  for (int k = 0; k < ngridMLSAtomList; k++)
606  {
607  dAdx(i,j) = dAdx(i,j) + dBdx(i,k) * P(k,j);
608  dAdy(i,j) = dAdy(i,j) + dBdy(i,k) * P(k,j);
609  dAdz(i,j) = dAdz(i,j) + dBdz(i,k) * P(k,j);
610  }
611  }
612  }
613 
614  for (int i = 0; i < 4; i++)
615  {
616  for (int j = 0; j < 4; j++)
617  {
618  for (int k = 0; k < 4; k++)
619  {
620  inverseAXdAdx(i,j) = inverseAXdAdx(i,j) + inverseA(i,k) * dAdx(k,j);
621  inverseAXdAdy(i,j) = inverseAXdAdy(i,j) + inverseA(i,k) * dAdy(k,j);
622  inverseAXdAdz(i,j) = inverseAXdAdz(i,j) + inverseA(i,k) * dAdz(k,j);
623  }
624  }
625  }
626 
627  for (int i = 0; i < 4; i++)
628  {
629  for (int j = 0; j < 4; j++)
630  {
631  for (int k = 0; k < 4; k++)
632  {
633  inverseAXdAdxXInverseA(i,j) = inverseAXdAdxXInverseA(i,j) + inverseAXdAdx(i,k) * inverseA(k,j);
634  inverseAXdAdyXInverseA(i,j) = inverseAXdAdyXInverseA(i,j) + inverseAXdAdy(i,k) * inverseA(k,j);
635  inverseAXdAdzXInverseA(i,j) = inverseAXdAdzXInverseA(i,j) + inverseAXdAdz(i,k) * inverseA(k,j);
636  }
637  }
638  }
639 
640  for (int i = 0; i < 4; i++)
641  {
642  for (int j = 0; j < ngridMLSAtomList; j++)
643  {
644  for (int k = 0; k < 4; k++)
645  {
646  fintmdx(i,j) = fintmdx(i,j) + inverseAXdAdxXInverseA(i,k) * B(k,j);
647  fintmdy(i,j) = fintmdy(i,j) + inverseAXdAdyXInverseA(i,k) * B(k,j);
648  fintmdz(i,j) = fintmdz(i,j) + inverseAXdAdzXInverseA(i,k) * B(k,j);
649  }
650  }
651  }
652 
653  for (int i = 0; i < 4; i++)
654  {
655  for (int j = 0; j < ngridMLSAtomList; j++)
656  {
657  for (int k = 0; k < 4; k++)
658  {
659  sintmdx(i,j) = sintmdx(i,j) + inverseA(i,k) * dBdx(k,j);
660  sintmdy(i,j) = sintmdy(i,j) + inverseA(i,k) * dBdy(k,j);
661  sintmdz(i,j) = sintmdz(i,j) + inverseA(i,k) * dBdz(k,j);
662  }
663  }
664  }
665 
666  for (int i = 0; i < 4; i++)
667  {
668  for (int j = 0; j < ngridMLSAtomList; j++)
669  {
670  fintmdx(i,j) = -fintmdx(i,j) + sintmdx(i,j);
671  fintmdy(i,j) = -fintmdy(i,j) + sintmdy(i,j);
672  fintmdz(i,j) = -fintmdz(i,j) + sintmdz(i,j);
673  }
674  }
675 
676  for (int i = 0; i < 4; i++)
677  {
678  for (int j = 0; j < 3; j++)
679  {
680  for (int k = 0; k < ngridMLSAtomList; k++)
681  {
682  dMlsDispdx(i,j) = dMlsDispdx(i,j) + fintmdx(i,k) * exactDisplacements(k,j);
683  dMlsDispdy(i,j) = dMlsDispdy(i,j) + fintmdy(i,k) * exactDisplacements(k,j);
684  dMlsDispdz(i,j) = dMlsDispdz(i,j) + fintmdz(i,k) * exactDisplacements(k,j);
685  }
686  }
687  }
688 
689  // debug
690 /*
691  std::cout << "fintmdx: " << std::endl;
692  std::cout << (fintmdx.transpose()).format(Eigen::FullPrecision) << std::endl;
693  std::cout << "fintmdy: " << std::endl;
694  std::cout << (fintmdy.transpose()).format(Eigen::FullPrecision) << std::endl;
695  std::cout << "fintmdz: " << std::endl;
696  std::cout << (fintmdz.transpose()).format(Eigen::FullPrecision) << std::endl;
697 
698  std::cout << "sintmdx: " << std::endl;
699  std::cout << (sintmdx.transpose()).format(Eigen::FullPrecision) << std::endl;
700  std::cout << "sintmdy: " << std::endl;
701  std::cout << (sintmdy.transpose()).format(Eigen::FullPrecision) << std::endl;
702  std::cout << "sintmdz: " << std::endl;
703  std::cout << (sintmdz.transpose()).format(Eigen::FullPrecision) << std::endl;
704 
705  std::cout << "dMlsDispdx: " << std::endl;
706  std::cout << (dMlsDispdx.transpose()).format(Eigen::FullPrecision) << std::endl;
707  std::cout << "dMlsDispdy: " << std::endl;
708  std::cout << (dMlsDispdy.transpose()).format(Eigen::FullPrecision) << std::endl;
709  std::cout << "dMlsDispdz: " << std::endl;
710  std::cout << (dMlsDispdz.transpose()).format(Eigen::FullPrecision) << std::endl;
711 */
712 
713  tensorF(0,0) = 1.0 + MlsDisp(1,0) \
714  + 1.0 * dMlsDispdx(0,0);
715  tensorF(1,0) = MlsDisp(2,0) \
716  + 1.0 * dMlsDispdy(0,0);
717  tensorF(2,0) = MlsDisp(3,0) \
718  + 1.0 * dMlsDispdz(0,0);
719  tensorF(0,1) = MlsDisp(1,1) \
720  + 1.0 * dMlsDispdx(0,1);
721  tensorF(1,1) = 1.0 + MlsDisp(2,1) \
722  + 1.0 * dMlsDispdy(0,1);
723  tensorF(2,1) = MlsDisp(3,1) \
724  + 1.0 * dMlsDispdz(0,1);
725  tensorF(0,2) = MlsDisp(1,2) \
726  + 1.0 * dMlsDispdx(0,2);
727  tensorF(1,2) = MlsDisp(2,2) \
728  + 1.0 * dMlsDispdy(0,2);
729  tensorF(2,2) = 1.0 + MlsDisp(3,2) \
730  + 1.0 * dMlsDispdz(0,2);
731 
732  gptIPushedF(0) = pgrid->coordinates[iGrid](0) + MlsDisp(0,0);
733  gptIPushedF(1) = pgrid->coordinates[iGrid](1) + MlsDisp(0,1);
734  gptIPushedF(2) = pgrid->coordinates[iGrid](2) + MlsDisp(0,2);
735  /*
736 
737  tensorF(0,0) = 1.0 + MlsDisp(1,0) \
738  + 1.0 * dMlsDispdx(0,0) + gridCoordinates[iGrid](0) * dMlsDispdx(1,0) \
739  + gridCoordinates[iGrid](1) * dMlsDispdx(2,0) + gridCoordinates[iGrid](2) * dMlsDispdx(3,0);
740  tensorF(0,1) = MlsDisp(2,0) \
741  + 1.0 * dMlsDispdy(0,0) + gridCoordinates[iGrid](0) * dMlsDispdy(1,0) \
742  + gridCoordinates[iGrid](1) * dMlsDispdy(2,0) + gridCoordinates[iGrid](2) * dMlsDispdy(3,0);
743  tensorF(0,2) = MlsDisp(3,0) \
744  + 1.0 * dMlsDispdz(0,0) + gridCoordinates[iGrid](0) * dMlsDispdz(1,0) \
745  + gridCoordinates[iGrid](1) * dMlsDispdz(2,0) + gridCoordinates[iGrid](2) * dMlsDispdz(3,0);
746  tensorF(1,0) = MlsDisp(1,1) \
747  + 1.0 * dMlsDispdx(0,1) + gridCoordinates[iGrid](0) * dMlsDispdx(1,1) \
748  + gridCoordinates[iGrid](1) * dMlsDispdx(2,1) + gridCoordinates[iGrid](2) * dMlsDispdx(3,1);
749  tensorF(1,1) = 1.0 + MlsDisp(2,1) \
750  + 1.0 * dMlsDispdy(0,1) + gridCoordinates[iGrid](0) * dMlsDispdy(1,1) \
751  + gridCoordinates[iGrid](1) * dMlsDispdy(2,1) + gridCoordinates[iGrid](2) * dMlsDispdy(3,1);
752  tensorF(1,2) = MlsDisp(3,1) \
753  + 1.0 * dMlsDispdz(0,1) + gridCoordinates[iGrid](0) * dMlsDispdz(1,1) \
754  + gridCoordinates[iGrid](1) * dMlsDispdz(2,1) + gridCoordinates[iGrid](2) * dMlsDispdz(3,1);
755  tensorF(2,0) = MlsDisp(1,2) \
756  + 1.0 * dMlsDispdx(0,2) + gridCoordinates[iGrid](0) * dMlsDispdx(1,2) \
757  + gridCoordinates[iGrid](1) * dMlsDispdx(2,2) + gridCoordinates[iGrid](2) * dMlsDispdx(3,2);
758  tensorF(2,1) = MlsDisp(2,2) \
759  + 1.0 * dMlsDispdy(0,2) + gridCoordinates[iGrid](0) * dMlsDispdy(1,2) \
760  + gridCoordinates[iGrid](1) * dMlsDispdy(2,2) + gridCoordinates[iGrid](2) * dMlsDispdy(3,2);
761  tensorF(2,2) = 1.0 + MlsDisp(3,2) \
762  + 1.0 * dMlsDispdz(0,2) + gridCoordinates[iGrid](0) * dMlsDispdz(1,2) \
763  + gridCoordinates[iGrid](1) * dMlsDispdz(2,2) + gridCoordinates[iGrid](2) * dMlsDispdz(3,2);
764 
765  gptIPushedF(0) = gridCoordinates[iGrid](0) \
766  + MlsDisp(0,0) + MlsDisp(1,0) * gridCoordinates[iGrid](0) \
767  + MlsDisp(2,0) * gridCoordinates[iGrid](1) \
768  + MlsDisp(3,0) * gridCoordinates[iGrid](2);
769  gptIPushedF(1) = gridCoordinates[iGrid](1) \
770  + MlsDisp(0,1) + MlsDisp(1,1) * gridCoordinates[iGrid](0) \
771  + MlsDisp(2,1) * gridCoordinates[iGrid](1) \
772  + MlsDisp(3,1) * gridCoordinates[iGrid](2);
773  gptIPushedF(2) = gridCoordinates[iGrid](2) \
774  + MlsDisp(0,2) + MlsDisp(1,2) * gridCoordinates[iGrid](0) \
775  + MlsDisp(2,2) * gridCoordinates[iGrid](1) \
776  + MlsDisp(3,2) * gridCoordinates[iGrid](2);
777 
778  */
779 
780  deformationGradient.push_back(tensorF);
781  gridPushed.push_back(gptIPushedF);
782  GridEnd:;
783 
784  //std::cout << "tensorF: " << std::endl;
785  //std::cout << tensorF.format(Eigen::FullPrecision) << std::endl;
786  //std::cout << "gptIPushedF: " << std::endl;
787  //std::cout << gptIPushedF.format(Eigen::FullPrecision) << std::endl;
788  }
789 }
790 
791 Mls::~Mls()
792 {
793  // TODO Auto-generated destructor stub
794 }
795 
796 void Mls::pushToCauchy(const std::vector<Matrix3d>& piolaStress,std::vector<Matrix3d>& cauchyStress)
797 {
798  MY_HEADING("Pushing the Piola-Kirchhoff Stress to the Cauchy stress.");
799  for (std::vector<Matrix3d>::size_type i = 0; i != piolaStress.size(); i++)
800  {
801  /*
802  // debug
803  Matrix3d A, B;
804  A(0,0) = 1;
805  A(1,0) = 2;
806  A(2,0) = 3;
807  A(0,1) = 4;
808  A(1,1) = 5;
809  A(2,1) = 6;
810  A(0,2) = 7;
811  A(1,2) = 8;
812  A(2,2) = 9;
813  std::cout << A << std::endl;
814 
815  B(0,0) = 9;
816  B(1,0) = 8;
817  B(2,0) = 7;
818  B(0,1) = 6;
819  B(1,1) = 5;
820  B(2,1) = 4;
821  B(0,2) = 3;
822  B(1,2) = 2;
823  B(2,2) = 1;
824  std::cout << B << std::endl;
825  std::cout << A * B.transpose() << std::endl;
826  // debug
827  */
828  cauchyStress.push_back(piolaStress[i] * deformationGradient[i].transpose() / deformationGradient[i].determinant());
829  }
830 }
831 
832 void Mls::writeDeformationGradient()
833 {
834  MY_HEADING("Writing Deformation Gradient.");
835  //std::setprecision(16);
836  std::ofstream file(name+".DeformationGradient");
837  file << deformationGradient.size() << "\n";
838  file << "\n";
839  for (auto& F : deformationGradient)
840  {
841 
842  Eigen::Map<Eigen::Matrix<double,1,DIM*DIM>> FRow(F.data(), F.size());
843  Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
844  file << FRow.format(fmt) << " " << std::setprecision(16) << F.determinant() << std::endl;
845  }
846 }
847 
848 void Mls::writeGridPushed()
849 {
850  MY_HEADING("Writing Pushed grid points.");
851  //std::setprecision(16);
852  std::ofstream file(name+".pushedGrids");
853  file << gridPushed.size() << "\n";
854  file << "\n";
855  for (auto& grid : gridPushed)
856  {
857  Eigen::Map<Eigen::Matrix<double,1,DIM>> gridRow(grid.data(), grid.size());
858  Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
859  file << gridRow.format(fmt) << std::endl;
860  }
861 }
862 
863 void Mls::writePushedCauchyStress(std::vector<Matrix3d>& cauchyStress)
864 {
865  MY_HEADING("Writing Pushed Cauchy Stress.");
866  //std::setprecision(16);
867  std::ofstream file(name+".CauchyPushed");
868  file << cauchyStress.size() << "\n";
869  file << "\n";
870  for (auto& pushedStress : cauchyStress)
871  {
872  Eigen::Map<Eigen::Matrix<double,1,DIM*DIM>> pushedStressRow(pushedStress.data(), pushedStress.size());
873  Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
874  file << pushedStressRow.format(fmt) << std::endl;
875  }
876 }
877 
878 
879 
880 
transpose
void transpose(double const *mat, double *const trans)
Definition: helper.hpp:59
neighbor_list.h
MY_HEADING
#define MY_HEADING(heading)
Definition: typedef.h:36
epsilon
const double epsilon
Definition: typedef.h:73
Mls::pushToCauchy
void pushToCauchy(const std::vector< Matrix3d > &piolaStress, std::vector< Matrix3d > &cauchyStress)
Definition: Mls.cpp:796
Configuration::coordinates
std::map< ConfigType, MatrixXd > coordinates
Definition: Configuration.h:24
MY_WARNING
#define MY_WARNING(message)
Definition: typedef.h:24
BoxConfiguration::getConfiguration
Configuration * getConfiguration(double) const
Definition: BoxConfiguration.cpp:90
Grid
Definition: Grid.h:31
Grid::getGridNeighborLists
std::vector< std::set< int > > getGridNeighborLists(const SubConfiguration &, const double &) const
Definition: Grid.cpp:83
VectorXd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor > VectorXd
Definition: typedef.h:59
GridBase::coordinates
std::vector< Vector3d > coordinates
Definition: Grid.h:18
Matrix4d
Eigen::Matrix< double, 4, 4, Eigen::RowMajor > Matrix4d
Definition: typedef.h:55
Mls::name
std::string name
Definition: Mls.h:18
MY_BANNER
#define MY_BANNER(announcement)
Definition: typedef.h:30
Stencil::expandStencil
void expandStencil(const Grid< configType > *pgrid, const double &contributingNeighborhoodSize, const double &noncontributingNeighborhoodSize)
Definition: Stencil.h:31
Mls::writePushedCauchyStress
void writePushedCauchyStress(std::vector< Matrix3d > &cauchyStress)
Definition: Mls.cpp:863
Mls::Mls
Mls(const BoxConfiguration &body, const Grid< Reference > *pgrid, double radiusMls, const std::string name)
Definition: Mls.cpp:21
MY_ERROR
#define MY_ERROR(message)
Definition: typedef.h:17
Vector3d
Eigen::Matrix< double, 1, DIM, Eigen::RowMajor > Vector3d
Definition: typedef.h:60
Mls::writeGridPushed
void writeGridPushed()
Definition: Mls.cpp:848
Matrix3d
Eigen::Matrix< double, DIM, DIM, Eigen::RowMajor > Matrix3d
Definition: typedef.h:56
Mls::gridPushed
std::vector< Vector3d > gridPushed
Definition: Mls.h:21
Mls::deformationGradient
std::vector< Matrix3d > deformationGradient
Definition: Mls.h:20
MatrixXd
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixXd
Definition: typedef.h:54
Mls::radiusMls
double radiusMls
Definition: Mls.h:19
Mls::writeDeformationGradient
void writeDeformationGradient()
Definition: Mls.cpp:832
Mls::~Mls
~Mls()
Definition: Mls.cpp:791