1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35
36 Christoph Kloss (DCS Computing GmbH, Linz)
37 Christoph Kloss (JKU Linz)
38 Philippe Seil (JKU Linz)
39
40 Copyright 2012- DCS Computing GmbH, Linz
41 Copyright 2009-2012 JKU Linz
42 ------------------------------------------------------------------------- */
43
44 /* ----------------------------------------------------------------------
45 Contributing author:
46 Evan Smuts (U Cape Town, surface velocity rotation)
47 ------------------------------------------------------------------------- */
48
49 #include "fix_mesh.h"
50 #include <stdio.h>
51 #include <string.h>
52 #include <algorithm>
53 #include <cmath>
54 #include "error.h"
55 #include "force.h"
56 #include "bounding_box.h"
57 #include "input_mesh_tri.h"
58 #include "fix_contact_history.h"
59 #include "fix_neighlist_mesh.h"
60 #include "tri_mesh_deform.h"
61 #include "fix_property_global.h"
62 #include "fix_move_mesh.h"
63 #include "tri_mesh_planar.h"
64 #include "modify.h"
65 #include "comm.h"
66 #include "math_extra.h"
67 #include "string_liggghts.h"
68
69 using namespace LAMMPS_NS;
70 using namespace FixConst;
71
72 #define EPSILON_V 0.00001
73
FixMesh(LAMMPS * lmp,int narg,char ** arg)74 FixMesh::FixMesh(LAMMPS *lmp, int narg, char **arg)
75 : FixBaseLiggghts(lmp, narg, arg),
76 atom_type_mesh_(-1),
77 temperature_mesh_(0.),
78 mass_temperature_(0.),
79 trackPerElementTemp_(false),
80 mesh_(NULL),
81 setupFlag_(false),
82 pOpFlag_(false),
83 manipulated_(false),
84 verbose_(false),
85 autoRemoveDuplicates_(false),
86 precision_(0.),
87 min_feature_length_(-1.),
88 element_exclusion_list_(0),
89 read_exclusion_list_(false),
90 exclusion_list_(0),
91 size_exclusion_list_(0),
92 fix_capacity_(0)
93 {
94 if(narg < 5)
95 error->fix_error(FLERR,this,"not enough arguments - at least keyword 'file' and a filename are required.");
96
97 do_support_multisphere();
98 do_not_need_radius();
99 do_not_need_mass();
100
101 restart_global = 1;
102
103 force_reneighbor = 1;
104 next_reneighbor = -1;
105
106 // parse args
107
108 iarg_ = 3;
109
110 char mesh_fname[256];
111 bool is_fix = strcmp(arg[iarg_], "fix") == 0;
112 if(!(strcmp(arg[iarg_],"file")==0 || is_fix))
113 error->fix_error(FLERR,this,"expecting keyword 'file' or 'fix'");
114 iarg_++;
115 strcpy(mesh_fname,arg[iarg_++]);
116
117 // parse args
118
119 bool hasargs = true;
120 while(iarg_ < narg && hasargs)
121 {
122 hasargs = false;
123 if(strcmp(arg[iarg_],"type") == 0) {
124 iarg_++;
125 atom_type_mesh_ = force->inumeric(FLERR,arg[iarg_++]);
126 if(atom_type_mesh_ < 1)
127 error->fix_error(FLERR,this,"'type' > 0 required");
128 hasargs = true;
129 } else if(strcmp(arg[iarg_],"verbose") == 0) {
130 if(narg < iarg_+2)
131 error->fix_error(FLERR,this,"not enough arguments for 'verbose'");
132 if(strcmp(arg[iarg_+1],"yes") == 0)
133 verbose_ = true;
134 else if(strcmp(arg[iarg_+1],"no"))
135 error->fix_error(FLERR,this,"expecing 'yes' or 'no' for 'verbose'");
136 iarg_ += 2;
137 hasargs = true;
138 } else if (strcmp(arg[iarg_],"region") == 0) {
139 if (iarg_+2 > narg) error->fix_error(FLERR,this,"not enough arguments for 'region'");
140 process_region(arg[iarg_+1]);
141 iarg_ += 2;
142 hasargs = true;
143 } else if(strcmp(arg[iarg_],"heal") == 0) {
144 if(narg < iarg_+2)
145 error->fix_error(FLERR,this,"not enough arguments for 'heal'");
146 if(strcmp(arg[iarg_+1],"auto_remove_duplicates") == 0)
147 autoRemoveDuplicates_ = true;
148 else if(strcmp(arg[iarg_+1],"no"))
149 error->fix_error(FLERR,this,"expecing 'auto_remove_duplicates' or 'no' for 'heal'");
150 iarg_ += 2;
151 hasargs = true;
152 } else if (strcmp(arg[iarg_],"precision") == 0) {
153 if (narg < iarg_+2) error->fix_error(FLERR,this,"not enough arguments");
154 iarg_++;
155 precision_ = force->numeric(FLERR,arg[iarg_++]);
156 if(precision_ < 0. || precision_ > 0.001)
157 error->fix_error(FLERR,this,"0 < precision < 0.001 required");
158 hasargs = true;
159 } else if (strcmp(arg[iarg_],"min_feature_length") == 0) {
160 if (narg < iarg_+2) error->fix_error(FLERR,this,"not enough arguments");
161 iarg_++;
162 min_feature_length_ = force->numeric(FLERR,arg[iarg_++]);
163 if(min_feature_length_ <= 0.)
164 error->fix_error(FLERR,this,"0 < min_feature_length > 0.0 required");
165 hasargs = true;
166 } else if (strcmp(arg[iarg_],"element_exclusion_list") == 0) {
167 if (narg < iarg_+3) error->fix_error(FLERR,this,"not enough arguments");
168 iarg_++;
169 if(strcmp("read",arg[iarg_]) == 0)
170 read_exclusion_list_ = true;
171 else if(strcmp("write",arg[iarg_]) == 0)
172 read_exclusion_list_ = false;
173 else error->fix_error(FLERR,this,"expecing 'read' or 'write' after 'element_exclusion_list'");
174 iarg_++;
175
176 if(1 < comm->nprocs && !read_exclusion_list_)
177 error->fix_error(FLERR,this,"'write' mode for 'element_exclusion_list' only works in serial");
178
179 // case write
180 if (!read_exclusion_list_ && 0 == comm->me)
181 {
182 element_exclusion_list_ = fopen(arg[iarg_++],"w");
183 if(!element_exclusion_list_)
184 error->one(FLERR,"Fix mesh: can not open file for 'element_exclusion_list' for writing");
185 }
186 // case read
187 else if(read_exclusion_list_ && 0 == comm->me)
188 {
189 element_exclusion_list_ = fopen(arg[iarg_++],"r");
190 if(!element_exclusion_list_)
191 error->one(FLERR,"Fix mesh: can not open file for 'element_exclusion_list' for reading");
192 }
193 if(0 < comm->me)
194 iarg_++;
195 hasargs = true;
196 }
197 }
198
199 if(min_feature_length_ > 0. && (!element_exclusion_list_ || read_exclusion_list_))
200 error->fix_error(FLERR,this,"'min_feature_length' requires use of 'element_exclusion_list write'");
201
202 // create/handle exclusion list
203 handle_exclusion_list();
204
205 // construct a mesh - can be surface or volume mesh
206 // just create object and return if reading data from restart file
207
208 if(modify->have_restart_data(this)) create_mesh_restart(mesh_fname);
209 else create_mesh(mesh_fname, is_fix);
210
211 // parse further args
212
213 hasargs = true;
214 while(iarg_ < narg && hasargs)
215 {
216 hasargs = false;
217 if(strcmp(arg[iarg_],"move") == 0) {
218 if (narg < iarg_+4) error->fix_error(FLERR,this,"not enough arguments");
219 moveMesh(force->numeric(FLERR,arg[iarg_+1]),force->numeric(FLERR,arg[iarg_+2]),force->numeric(FLERR,arg[iarg_+3]));
220 manipulated_ = true;
221 iarg_ += 4;
222 hasargs = true;
223 } else if(strcmp(arg[iarg_],"rotate") == 0) {
224 if (narg < iarg_+7)
225 error->fix_error(FLERR,this,"not enough arguments");
226 if(strcmp(arg[iarg_+1],"axis"))
227 error->fix_error(FLERR,this,"expecting keyword 'axis' after keyword 'rotate'");
228 if(strcmp(arg[iarg_+5],"angle"))
229 error->fix_error(FLERR,this,"expecting keyword 'angle' after axis definition");
230 rotateMesh(force->numeric(FLERR,arg[iarg_+2]),force->numeric(FLERR,arg[iarg_+3]),force->numeric(FLERR,arg[iarg_+4]),
231 force->numeric(FLERR,arg[iarg_+6]));
232 manipulated_ = true;
233 iarg_ += 7;
234 hasargs = true;
235 } else if(strcmp(arg[iarg_],"scale") == 0) {
236 if (narg < iarg_+2) error->fix_error(FLERR,this,"not enough arguments");
237 scaleMesh(force->numeric(FLERR,arg[iarg_+1]));
238 manipulated_ = true;
239 iarg_ += 2;
240 hasargs = true;
241 } else if (strcmp(arg[iarg_],"temperature") == 0) {
242 iarg_++;
243 temperature_mesh_ = atof(arg[iarg_++]);
244 mesh_->prop().addGlobalProperty< ScalarContainer<double> >("Temp","comm_none","frame_invariant","restart_yes");
245 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("Temp",temperature_mesh_);
246 mesh_->prop().addGlobalProperty< ScalarContainer<double> >("heatFlux","comm_none","frame_invariant","restart_no");
247 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("heatFlux",0.);
248 mesh_->prop().addGlobalProperty< ScalarContainer<double> >("heatFluxTotal","comm_none","frame_invariant","restart_yes");
249 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("heatFluxTotal",0.);
250
251 hasargs = true;
252 } else if (strcmp(arg[iarg_],"mass_temperature") == 0) {
253 iarg_++;
254 mass_temperature_ = atof(arg[iarg_++]);
255 if(mass_temperature_ <= 0.)
256 error->fix_error(FLERR,this,"mass_temperature > 0 expected");
257 hasargs = true;
258 }
259 }
260 }
261
262 /* ---------------------------------------------------------------------- */
263
handle_exclusion_list()264 void FixMesh::handle_exclusion_list()
265 {
266 // case read exclusion list
267 if(read_exclusion_list_)
268 {
269 // read from exclusion list, only on proc 0
270 if(element_exclusion_list_)
271 {
272 char read_string[200];
273 while(fgets(read_string,200,element_exclusion_list_ ) != 0)
274 {
275 // remove trailing newline
276 read_string[strcspn(read_string, "\r\n")] = 0; // works for LF, CR, CRLF, LFCR, ...
277 // trim string
278 strtrim(read_string);
279
280 int line = force->inumeric(FLERR,read_string);
281 memory->grow(exclusion_list_, size_exclusion_list_+1, "exclusion_list");
282 exclusion_list_[size_exclusion_list_++] = line;
283 }
284 }
285 // send size_exclusion_list_ to all procs
286 MPI_Max_Scalar(size_exclusion_list_,world);
287 if(0 < comm->me)
288 {
289 memory->grow(exclusion_list_, size_exclusion_list_, "exclusion_list");
290 vectorZeroizeN(exclusion_list_,size_exclusion_list_);
291 }
292 MPI_Max_Vector(exclusion_list_,size_exclusion_list_,world);
293
294 // sort
295 if(size_exclusion_list_ > 0)
296 {
297 std::vector<int> sorted;
298 for(int i = 0; i < size_exclusion_list_;i++)
299 sorted.push_back(exclusion_list_[i]);
300 sort(sorted.begin(),sorted.end());
301 for(int i = 0; i < size_exclusion_list_;i++)
302 exclusion_list_[i] = sorted[i];
303 }
304 }
305 }
306
307 /* ---------------------------------------------------------------------- */
308
~FixMesh()309 FixMesh::~FixMesh()
310 {
311 delete mesh_;
312 if (element_exclusion_list_)
313 fclose(element_exclusion_list_);
314
315 if(exclusion_list_)
316 memory->sfree(exclusion_list_);
317 }
318
319 /* ---------------------------------------------------------------------- */
320
post_create()321 void FixMesh::post_create()
322 {
323 // check if all element property container have same length
324 // could potentially be whacked by adding element properties
325 // at the wrong place in code
326 mesh_->check_element_property_consistency();
327
328 }
329
330 /* ---------------------------------------------------------------------- */
331
create_mesh(char * mesh_fname,bool is_fix)332 void FixMesh::create_mesh(char *mesh_fname, bool is_fix)
333 {
334
335 if(strncmp(style,"mesh/surface",12) == 0)
336 {
337 if(strcmp(style,"mesh/surface/stress/deform") == 0)
338 mesh_ = new TriMeshDeformable(lmp);
339 else if(strcmp(style,"mesh/surface/planar") == 0)
340 mesh_ = new TriMeshPlanar(lmp);
341 else
342 mesh_ = new TriMesh(lmp);
343
344 // set properties that are important for reading
345 mesh_->setMeshID(id);
346 if(verbose_) mesh_->setVerbose();
347 if(autoRemoveDuplicates_) mesh_->autoRemoveDuplicates();
348 if(precision_ > 0.) mesh_->setPrecision(precision_);
349 if(min_feature_length_ > 0.) mesh_->setMinFeatureLength(min_feature_length_);
350
351 if (is_fix)
352 {
353 Fix *fix_base = modify->find_fix_id(mesh_fname);
354 if (!fix_base)
355 error->all(FLERR, "Could not find appropriate fix to read mesh data from");
356 if (!fix_base->can_create_mesh())
357 error->all(FLERR, "Supplied fix does not provide a mesh");
358 const int tri_count = fix_base->getCreateMeshTriCount();
359 if (tri_count == 0)
360 error->all(FLERR, "Extruded triangle count is zero. Are you sure your fix mesh/surface has extrude_planar set");
361 double **node;
362 node = new double*[3];
363 for (int i=0; i < tri_count*3; i++)
364 {
365 node[i%3] = fix_base->getCreateMeshTriNode(i);
366 if (i%3 == 2)
367 static_cast<TriMesh*>(mesh_)->addElement(node, -1);
368 }
369 delete [] node;
370 }
371 else
372 {
373 // read file
374 // can be from STL file or VTK file
375 InputMeshTri *mesh_input = new InputMeshTri(lmp,0,NULL);
376
377 // case write exlusion list
378 if(!read_exclusion_list_ && element_exclusion_list_)
379 mesh_->setElementExclusionList(element_exclusion_list_);
380
381 mesh_input->meshtrifile(mesh_fname,static_cast<TriMesh*>(mesh_),verbose_,
382 size_exclusion_list_,exclusion_list_,region_);
383
384 delete mesh_input;
385 }
386 }
387 else error->one(FLERR,"Illegal implementation of create_mesh();");
388 }
389
390 /* ---------------------------------------------------------------------- */
391
create_mesh_restart(char * mesh_fname)392 void FixMesh::create_mesh_restart(char *mesh_fname)
393 {
394
395 if(strcmp(style,"mesh/surface/stress/deform") == 0)
396 mesh_ = new TriMeshDeformable(lmp);
397 else if(strcmp(style,"mesh/surface/planar") == 0)
398 mesh_ = new TriMeshPlanar(lmp);
399 else if(strncmp(style,"mesh/surface",12) == 0)
400 mesh_ = new TriMesh(lmp);
401 else error->one(FLERR,"Illegal implementation of create_mesh();");
402
403 // test if file exists
404 if (0 == comm->me)
405 {
406 char str[512];
407 FILE * test = fopen(mesh_fname,"r");
408 if(!test)
409 {
410
411 sprintf(str,"Cannot open mesh file %s. FYI: This file is required, but data will be taken from restart file",mesh_fname);
412 error->one(FLERR,str);
413 }
414 else
415 {
416 sprintf(str,"INFO: mesh file (%s) is required, but data will be taken from restart file",mesh_fname);
417 error->message(FLERR,str);
418 fclose(test);
419 }
420 }
421
422 // set properties that are important for reading
423 mesh_->setMeshID(id);
424 if(verbose_) mesh_->setVerbose();
425 if(autoRemoveDuplicates_) mesh_->autoRemoveDuplicates();
426 if(precision_ > 0.) mesh_->setPrecision(precision_);
427
428 // case write exlusion list
429 if(!read_exclusion_list_ && element_exclusion_list_)
430 mesh_->setElementExclusionList(element_exclusion_list_);
431 }
432
433 /* ---------------------------------------------------------------------- */
434
pre_delete(bool unfixflag)435 void FixMesh::pre_delete(bool unfixflag)
436 {
437 // error if moving mesh is operating on a mesh to be deleted
438
439 // also error if dump is operating on mesh
440 if(unfixflag)
441 {
442 const int n_move = modify->n_fixes_style("move/mesh");
443 if(mesh_->isMoving() && n_move > 0)
444 {
445 for (int i = 0; i < n_move; i++) {
446 const FixMoveMesh* fix_move = static_cast<FixMoveMesh*>(modify->find_fix_style_strict("move/mesh", i));
447 if (fix_move->fixMeshCompare(this))
448 {
449 std::string error_msg =
450 std::string("illegal unfix command, may not unfix a moving mesh while a fix move is applied to it. ") +
451 "Unfix the fix move/mesh first (id: " +
452 fix_move->id +
453 ")";
454 error->fix_error(FLERR,this,error_msg.c_str());
455 }
456 }
457 }
458 }
459 }
460
461 /* ---------------------------------------------------------------------- */
462
init()463 void FixMesh::init()
464 {
465 FixBaseLiggghts::init();
466
467 if(mass_temperature_ > 0.)
468 {
469 int max_type = atom->get_properties()->max_type();
470 fix_capacity_ =
471 static_cast<FixPropertyGlobal*>(modify->find_fix_property("thermalCapacity","property/global","peratomtype",max_type,0,style));
472
473 }
474 }
475
476 /* ---------------------------------------------------------------------- */
477
setmask()478 int FixMesh::setmask()
479 {
480 int mask = 0;
481 mask |= PRE_EXCHANGE;
482 mask |= PRE_FORCE;
483 mask |= FINAL_INTEGRATE;
484 return mask;
485 }
486
487 /* ---------------------------------------------------------------------- */
488
setup_pre_force(int vflag)489 void FixMesh::setup_pre_force(int vflag)
490 {
491 // first-time set-up
492
493 if(!setupFlag_)
494 {
495 initialSetup();
496 setupFlag_ = true;
497 }
498 // if mesh already set-up and parallelized
499
500 else
501 {
502 mesh_->pbcExchangeBorders(1);
503 }
504
505 // clear reverse comm properties
506 mesh_->clearReverse();
507
508 pOpFlag_ = false;
509
510 }
511
512 /* ---------------------------------------------------------------------- */
513
setup(int vflag)514 void FixMesh::setup(int vflag)
515 {
516 if(temperature_mesh_ > 1e-12)
517 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("Temp",temperature_mesh_);
518
519 mesh_->reverseComm();
520 }
521
522 /* ---------------------------------------------------------------------- */
523
initialSetup()524 void FixMesh::initialSetup()
525 {
526
527 mesh_->initialSetup();
528
529 // warn if there are elements that extend outside box
530 if(!mesh_->allNodesInsideSimulationBox() && 0 == comm->me)
531 error->warning(FLERR,"Not all nodes of fix mesh inside simulation box, "
532 "elements will be deleted or wrapped around periodic boundary conditions");
533 if(comm->me == 0)
534 fprintf(screen,"Import and parallelization of mesh %s containing %d triangle(s) successful\n",
535 id,mesh_->sizeGlobal());
536 }
537
538 /* ----------------------------------------------------------------------
539 invoke parallelism
540 ------------------------------------------------------------------------- */
541
pre_exchange()542 void FixMesh::pre_exchange()
543 {
544
545 pOpFlag_ = true;
546 }
547
548 /* ----------------------------------------------------------------------
549 forward comm for mesh
550 ------------------------------------------------------------------------- */
551
pre_force(int vflag)552 void FixMesh::pre_force(int vflag)
553 {
554 // case re-neigh step
555 if(pOpFlag_)
556 {
557
558 mesh_->pbcExchangeBorders(0);
559
560 pOpFlag_ = false;
561 }
562 // case regular step
563 else
564 {
565 mesh_->forwardComm();
566
567 if(mesh_->decideRebuild())
568 {
569
570 next_reneighbor = update->ntimestep + 1;
571 }
572 }
573
574 // clear reverse comm properties
575 mesh_->clearReverse();
576 }
577
578 /* ----------------------------------------------------------------------
579 reverse comm for mesh
580 ------------------------------------------------------------------------- */
581
final_integrate()582 void FixMesh::final_integrate()
583 {
584
585 mesh_->reverseComm();
586
587 bool has_per_element_heattransfer = (0 == strcmp("mesh/surface",style) && 0 == strcmp("heattransfer", style));
588
589 if(!has_per_element_heattransfer && mass_temperature_ > 0. && mesh_->prop().getGlobalProperty< ScalarContainer<double> >("Temp"))
590 {
591 double Temp_wall = (*mesh_->prop().getGlobalProperty< ScalarContainer<double> >("Temp"))(0);
592 double flux = (*mesh_->prop().getGlobalProperty< ScalarContainer<double> >("heatFlux"))(0);
593 MPI_Sum_Scalar(flux,world);
594 double dt = update->dt;
595
596 double capacity = fix_capacity_->compute_vector(atom_type_mesh_-1);
597 Temp_wall += flux * dt / (mass_temperature_*capacity);
598
599 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("Temp",Temp_wall);
600 mesh_->prop().setGlobalProperty< ScalarContainer<double> >("heatFlux",0.);
601 }
602 }
603
604 /* ---------------------------------------------------------------------- */
605
min_type()606 int FixMesh::min_type()
607 {
608 if(atom_type_mesh_ == -1) return 1;
609 return atom_type_mesh_;
610 }
611
612 /* ---------------------------------------------------------------------- */
613
max_type()614 int FixMesh::max_type()
615 {
616 if(atom_type_mesh_ == -1) return 1;
617 return atom_type_mesh_;
618 }
619
620 /* ----------------------------------------------------------------------
621 moves the mesh by a vector - (dx dy dz) is the displacement vector
622 ------------------------------------------------------------------------- */
623
moveMesh(double const dx,double const dy,double const dz)624 void FixMesh::moveMesh(double const dx, double const dy, double const dz)
625 {
626 if (comm->me == 0)
627 {
628 //fprintf(screen,"moving mesh by ");
629 //fprintf(screen,"%f %f %f\n", dx,dy,dz);
630 }
631
632 const double arg[3] = {dx,dy,dz};
633 mesh_->move(arg);
634 }
635
636 /* ----------------------------------------------------------------------
637 rotates the mesh around an axis through the origin
638 phi the rotation angle
639 (axisX axisY axisZ) vector in direction of the axis
640 ------------------------------------------------------------------------- */
641
rotateMesh(double const axisX,double const axisY,double const axisZ,double const phi)642 void FixMesh::rotateMesh(double const axisX, double const axisY, double const axisZ, double const phi)
643 {
644 double axis[3] = {axisX,axisY,axisZ}, p[3] = {0.,0.,0.};
645
646 if(vectorMag3D(axis) < 1e-5)
647 error->fix_error(FLERR,this,"illegal magnitude of rotation axis");
648
649 if (comm->me == 0)
650 {
651 //fprintf(screen,"rotate ");
652 //fprintf(screen,"%f %f %f %f\n", phi, axisX, axisY, axisZ);
653 }
654
655 mesh_->rotate(phi*3.14159265/180.0,axis,p);
656 }
657
658 /* ----------------------------------------------------------------------
659 scales the mesh by a factor in x, y and z direction
660 can also be used to distort a mesh
661 (factorX factorY factorZ) vector contains the factors to scale the
662 mesh in the xyz direction
663 ------------------------------------------------------------------------- */
664
scaleMesh(double const factor)665 void FixMesh::scaleMesh(double const factor)
666 {
667 if (comm->me == 0){
668 //fprintf(screen,"scale ");
669 //fprintf(screen,"%f \n", factor);
670 }
671 mesh_->scale(factor);
672 }
673
674 /* ----------------------------------------------------------------------
675 pack entire state of Fix into one write
676 ------------------------------------------------------------------------- */
677
write_restart(FILE * fp)678 void FixMesh::write_restart(FILE *fp)
679 {
680 mesh_->writeRestart(fp);
681 }
682
683 /* ----------------------------------------------------------------------
684 use state info from restart file to restart the Fix
685 ------------------------------------------------------------------------- */
686
restart(char * buf)687 void FixMesh::restart(char *buf)
688 {
689 double *list = (double *) buf;
690 mesh_->restart(list);
691 }
692
693 /* ----------------------------------------------------------------------
694 change box extent due to mesh node position
695 ------------------------------------------------------------------------- */
696
box_extent(double & xlo,double & xhi,double & ylo,double & yhi,double & zlo,double & zhi)697 void FixMesh::box_extent(double &xlo,double &xhi,double &ylo,double &yhi,double &zlo,double &zhi)
698 {
699 double node[3];
700 int size = mesh()->size();
701 int numNodes = mesh()->numNodes();
702
703 for(int i = 0; i < size; i++)
704 {
705
706 for(int j = 0; j < numNodes; j++)
707 {
708 mesh()->node_slow(i,j,node);
709 xlo = std::min(xlo,node[0]);
710 xhi = std::max(xhi,node[0]);
711 ylo = std::min(ylo,node[1]);
712 yhi = std::max(yhi,node[1]);
713 zlo = std::min(zlo,node[2]);
714 zhi = std::max(zhi,node[2]);
715 }
716 }
717 }
718
719 /* ---------------------------------------------------------------------- */
720
move(const double * const dx,const FixMoveMesh * const caller)721 void FixMesh::move(const double * const dx, const FixMoveMesh * const caller)
722 {
723 mesh_->move(dx);
724 std::list<FixMoveMesh *>::iterator it;
725 bool found = false;
726 for (it = fixMoveMeshes_.begin(); it != fixMoveMeshes_.end(); it++)
727 {
728 if (*it == caller)
729 found = true;
730 if (found)
731 (*it)->move(dx);
732 }
733 }
734
735 /* ---------------------------------------------------------------------- */
736
rotate(const double dphi,const double * const axis,const double * const center,const FixMoveMesh * const caller)737 void FixMesh::rotate(const double dphi, const double * const axis, const double * const center, const FixMoveMesh * const caller)
738 {
739 mesh_->rotate(dphi, axis, center);
740 std::list<FixMoveMesh *>::iterator it;
741 bool found = false;
742 for (it = fixMoveMeshes_.begin(); it != fixMoveMeshes_.end(); it++)
743 {
744 if (*it == caller)
745 found = true;
746 if (found)
747 (*it)->rotate(dphi, axis, center);
748 }
749 }
750