MDStressLab++
Loading...
Searching...
No Matches
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)
21Mls::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 GridSubConfiguration<Reference> neighborListsOfReferenceGridsMls(*pgrid,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 = neighborListsOfReferenceGridsMls.getGridPointNeighbors(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
791Mls::~Mls()
792{
793 // TODO Auto-generated destructor stub
794}
795
796void 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
832void 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
840 Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "");
841 file << std::fixed << std::setprecision(std::numeric_limits<double>::max_digits10);
842 file << "Properties=pos:R:3:DeformationGradient:R:9" << std::endl;
843 int index= 0;
844 //XX, YX, ZX, XY, YY, ZY, XZ, YZ, ZZ
845 for (auto& F : deformationGradient)
846 {
847 file << gridPushed[index]
848 << std::setw(25) << F(0,0)
849 << std::setw(25) << F(1,0)
850 << std::setw(25) << F(2,0)
851 << std::setw(25) << F(0,1)
852 << std::setw(25) << F(1,1)
853 << std::setw(25) << F(2,1)
854 << std::setw(25) << F(0,2)
855 << std::setw(25) << F(1,2)
856 << std::setw(25) << F(2,2)
857 << std::endl;
858 index++;
859
860 }
861 /*
862 for (auto& F : deformationGradient)
863 {
864
865 Eigen::Map<Eigen::Matrix<double,1,DIM*DIM>> FRow(F.data(), F.size());
866 Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
867 file << FRow.format(fmt) << " " << std::setprecision(16) << F.determinant() << std::endl;
868 }
869 */
870}
871
872void Mls::writeGridPushed()
873{
874 MY_HEADING("Writing Pushed grid points.");
875 //std::setprecision(16);
876 std::ofstream file(name+"Pushed.grid");
877 file << gridPushed.size() << "\n";
878 file << "\n";
879 for (auto& grid : gridPushed)
880 {
881 Eigen::Map<Eigen::Matrix<double,1,DIM>> gridRow(grid.data(), grid.size());
882 Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
883 file << gridRow.format(fmt) << std::endl;
884 }
885}
886
887void Mls::writePushedCauchyStress(std::vector<Matrix3d>& cauchyStress)
888{
889 MY_HEADING("Writing Pushed Cauchy Stress.");
890 /*
891 std::ofstream file(name+".CauchyPushed");
892 file << cauchyStress.size() << "\n";
893 for (auto& pushedStress : cauchyStress)
894 {
895 Eigen::Map<Eigen::Matrix<double,1,DIM*DIM>> pushedStressRow(pushedStress.data(), pushedStress.size());
896 Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "", "");
897 file << pushedStressRow.format(fmt) << std::endl;
898 }
899 */
900
901 std::ofstream file(name+"Pushed.stress");
902 file << cauchyStress.size() << "\n";
903 Eigen::IOFormat fmt(Eigen::FullPrecision, 0, " ", "\n", "", "", "");
904 file << std::fixed << std::setprecision(std::numeric_limits<double>::max_digits10);
905 file << "Properties=pos:R:3:stress:R:6" << std::endl;
906 int index= 0;
907 for (auto& pushedStress : cauchyStress)
908 {
909 file << gridPushed[index]
910 << std::setw(25) << pushedStress(0,0)
911 << std::setw(25) << pushedStress(1,1)
912 << std::setw(25) << pushedStress(2,2)
913 << std::setw(25) << pushedStress(0,1)
914 << std::setw(25) << pushedStress(0,2)
915 << std::setw(25) << pushedStress(1,2)
916 << std::endl;
917 index++;
918
919 }
920}
921
922
923
924
BoxConfiguration
Represents a particle configuration including simulation box information.
Definition BoxConfiguration.h:33
BoxConfiguration::pbc
Vector3i pbc
Periodic boundary conditions. pbc=(1,0,1) implies periodicity along the and -directions.
Definition BoxConfiguration.h:43
BoxConfiguration::getConfiguration
Configuration * getConfiguration(double padding) const
This function returns a padded configuration by adding padding atoms originating due to pbcs.
Definition BoxConfiguration.cpp:275
BoxConfiguration::box
Matrix3d box
The current and reference box vectors stored as columns of respective matrices.
Definition BoxConfiguration.h:38
GridBase::coordinates
std::vector< Vector3d > coordinates
Definition Grid.h:18
GridSubConfiguration::getGridPointNeighbors
std::set< int > getGridPointNeighbors(const int &) const
Definition Grid.cpp:174
Grid
Definition Grid.h:31
Mls::Mls
Mls(const BoxConfiguration &body, const Grid< Reference > *pgrid, double radiusMls, const std::string name)
Constructs the deformation gradient field from a body configuration and grid.
Definition Mls.cpp:21
Mls::~Mls
~Mls()
Definition Mls.cpp:791
Mls::deformationGradient
std::vector< Matrix3d > deformationGradient
Deformation gradient at each grid point (one per grid coordinate).
Definition Mls.h:42
Mls::writeDeformationGradient
void writeDeformationGradient()
Writes the deformation gradient at each grid point to a file.
Definition Mls.cpp:832
Mls::pushToCauchy
void pushToCauchy(const std::vector< Matrix3d > &piolaStress, std::vector< Matrix3d > &cauchyStress)
Converts Piola-Kirchhoff stress to Cauchy stress at grid points.
Definition Mls.cpp:796
Mls::radiusMls
double radiusMls
Radius of influence for the MLS weighting function.
Definition Mls.h:37
Mls::name
std::string name
Identifier name used for output file prefixes.
Definition Mls.h:32
Mls::writePushedCauchyStress
void writePushedCauchyStress(std::vector< Matrix3d > &cauchyStress)
Writes Cauchy stresses at grid points to a file.
Definition Mls.cpp:887
Mls::gridPushed
std::vector< Vector3d > gridPushed
Positions of grid points after applying MLS displacements.
Definition Mls.h:47
Mls::writeGridPushed
void writeGridPushed()
Writes the pushed (deformed) grid point coordinates to a file.
Definition Mls.cpp:872
Stencil::expandStencil
void expandStencil(const std::vector< Vector3d > &gridCoordinates, const MatrixXd &coordinates, const double &contributingNeighborhoodSize, const double &noncontributingNeighborhoodSize)
Definition Stencil.h:30
transpose
void transpose(double const *mat, double *const trans)
Definition helper.hpp:59
neighbor_list.h
MY_ERROR
#define MY_ERROR(message)
Definition typedef.h:17
Matrix3d
Eigen::Matrix< double, DIM, DIM, Eigen::RowMajor > Matrix3d
Definition typedef.h:56
MatrixXd
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixXd
Definition typedef.h:54
VectorXd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor > VectorXd
Definition typedef.h:59
Vector3d
Eigen::Matrix< double, 1, DIM, Eigen::RowMajor > Vector3d
Definition typedef.h:60
Matrix4d
Eigen::Matrix< double, 4, 4, Eigen::RowMajor > Matrix4d
Definition typedef.h:55
Current
@ Current
Definition typedef.h:70
Reference
@ Reference
Definition typedef.h:69
epsilon
const double epsilon
Definition typedef.h:73
MY_HEADING
#define MY_HEADING(heading)
Definition typedef.h:36
MY_WARNING
#define MY_WARNING(message)
Definition typedef.h:24
MY_BANNER
#define MY_BANNER(announcement)
Definition typedef.h:30