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 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Christoph Kloss (DCS Computing GmbH, Linz)
39 Christoph Kloss (JKU Linz)
40 Philippe Seil (JKU Linz)
41 Arno Mayrhofer (CFDEMresearch GmbH, Linz)
42
43 Copyright 2012- DCS Computing GmbH, Linz
44 Copyright 2009-2012 JKU Linz
45 Copyright 2016- CFDEMresearch GmbH, Linz
46 ------------------------------------------------------------------------- */
47
48 /* ----------------------------------------------------------------------
49 Contributing author, surface velocity rotation only:
50 Evan Smuts (U Cape Town)
51 ------------------------------------------------------------------------- */
52
53 #include "fix_mesh_surface.h"
54 #include <stdio.h>
55 #include <string>
56 #include <sstream>
57 #include "error.h"
58 #include "group.h"
59 #include "force.h"
60 #include "bounding_box.h"
61 #include "input_mesh_tri.h"
62 #include "fix_contact_history_mesh.h"
63 #include "fix_contact_property_atom_wall.h"
64 #include "fix_neighlist_mesh.h"
65 #include "multi_node_mesh.h"
66 #include "modify.h"
67 #include "comm.h"
68 #include "math_extra.h"
69 #include "style_mesh_module.h"
70 #include "mesh_module_stress.h"
71 #include "variable.h"
72
73 using namespace LAMMPS_NS;
74 using namespace FixConst;
75
76 #define EPSILON_V 0.00001
77
78 enum{NONE,CONSTANT,EQUAL};
79
FixMeshSurface(LAMMPS * lmp,int narg,char ** arg)80 FixMeshSurface::FixMeshSurface(LAMMPS *lmp, int narg, char **arg)
81 : FixMesh(lmp, narg, arg),
82 fix_contact_history_mesh_(0),
83 fix_mesh_neighlist_(0),
84 fix_meshforce_contact_(NULL),
85 fix_meshforce_contact_stress_(NULL),
86 fix_mesh_multicontact_data_(NULL),
87 velFlag_(false),
88 vSurfStrX_(NULL),
89 vSurfStrY_(NULL),
90 vSurfStrZ_(NULL),
91 vSurfVarX_(-1),
92 vSurfVarY_(-1),
93 vSurfVarZ_(-1),
94 vSurfStyleX_(-1),
95 vSurfStyleY_(-1),
96 vSurfStyleZ_(-1),
97 angVelFlag_(false),
98 omegaStr_(NULL),
99 omegaVar_(-1),
100 omegaStyle_(-1),
101 n_dump_active_(0),
102 curvature_(0.),
103 curvature_tolerant_(false),
104 extrude_mesh_(false),
105 extrusion_length_(0.0),
106 extrusion_tri_count_(0),
107 extrusion_tri_nodes_(NULL),
108 extrusion_created_(false)
109 {
110 // check if type has been read
111 if(atom_type_mesh_ == -1)
112 error->fix_error(FLERR,this,"expecting keyword 'type'");
113
114 // parse further args
115
116 bool hasargs = true;
117 while(iarg_ < narg && hasargs)
118 {
119 hasargs = false;
120 if (strcmp(arg[iarg_],"surface_vel") == 0) {
121 if (narg < iarg_+4) error->fix_error(FLERR,this,"not enough arguments");
122 iarg_++;
123 velFlag_ = true;
124 if (strstr(arg[iarg_], "v_") == arg[iarg_])
125 {
126 int n = strlen(&arg[iarg_][2]) + 1;
127 vSurfStrX_ = new char[n];
128 strcpy(vSurfStrX_, &arg[iarg_++][2]);
129 }
130 else
131 {
132 vSurf_[0] = force->numeric(FLERR,arg[iarg_++]);
133 vSurfStyleX_ = CONSTANT;
134 }
135 if (strstr(arg[iarg_], "v_") == arg[iarg_])
136 {
137 int n = strlen(&arg[iarg_][2]) + 1;
138 vSurfStrY_ = new char[n];
139 strcpy(vSurfStrY_, &arg[iarg_++][2]);
140 }
141 else
142 {
143 vSurf_[1] = force->numeric(FLERR,arg[iarg_++]);
144 vSurfStyleY_ = CONSTANT;
145 }
146 if (strstr(arg[iarg_], "v_") == arg[iarg_])
147 {
148 int n = strlen(&arg[iarg_][2]) + 1;
149 vSurfStrZ_ = new char[n];
150 strcpy(vSurfStrZ_, &arg[iarg_++][2]);
151 }
152 else
153 {
154 vSurf_[2] = force->numeric(FLERR,arg[iarg_++]);
155 vSurfStyleZ_ = CONSTANT;
156 }
157 hasargs = true;
158 } else if (strcmp(arg[iarg_],"surface_ang_vel") == 0) {
159 if (narg < iarg_+11) error->fix_error(FLERR,this,"not enough arguments");
160 iarg_++;
161 angVelFlag_ = true;
162 if(strcmp(arg[iarg_++],"origin"))
163 error->fix_error(FLERR,this,"expecting keyword 'origin' after 'rotation'");
164 origin_[0] = force->numeric(FLERR,arg[iarg_++]);
165 origin_[1] = force->numeric(FLERR,arg[iarg_++]);
166 origin_[2] = force->numeric(FLERR,arg[iarg_++]);
167 if(strcmp(arg[iarg_++],"axis"))
168 error->fix_error(FLERR,this,"expecting keyword 'axis' after definition of 'origin'");
169 axis_[0] = force->numeric(FLERR,arg[iarg_++]);
170 axis_[1] = force->numeric(FLERR,arg[iarg_++]);
171 axis_[2] = force->numeric(FLERR,arg[iarg_++]);
172 if(vectorMag3D(axis_) < EPSILON_V)
173 error->fix_error(FLERR,this,"'axis' vector must me non-zero");
174 if(strcmp(arg[iarg_++],"omega"))
175 error->fix_error(FLERR,this,"expecting keyword 'omega' after definition of 'axis'");
176 // positive omega give anti-clockwise (CCW) rotation
177 if (strstr(arg[iarg_], "v_") == arg[iarg_])
178 {
179 int n = strlen(&arg[iarg_][2]) + 1;
180 omegaStr_ = new char[n];
181 strcpy(omegaStr_, &arg[iarg_++][2]);
182 }
183 else
184 {
185 omegaSurf_ = force->numeric(FLERR,arg[iarg_++]);
186 omegaStyle_ = CONSTANT;
187 }
188 hasargs = true;
189 } else if (strcmp(arg[iarg_],"curvature") == 0) {
190 if (narg < iarg_+2)
191 error->fix_error(FLERR,this,"not enough arguments for 'curvature'");
192 iarg_++;
193 curvature_ = force->numeric(FLERR,arg[iarg_++]);
194 if(curvature_ <= 0. || curvature_ > 60)
195 error->fix_error(FLERR,this,"0° < curvature < 60° required");
196 curvature_ = cos(curvature_*M_PI/180.);
197 hasargs = true;
198 } else if (strcmp(arg[iarg_],"curvature_tolerant") == 0) {
199 if (narg < iarg_+2)
200 error->fix_error(FLERR,this,"not enough arguments for 'curvature_tolerant'");
201 iarg_++;
202 if(0 == strcmp(arg[iarg_],"yes"))
203 curvature_tolerant_= true;
204 else if(0 == strcmp(arg[iarg_],"no"))
205 curvature_tolerant_= false;
206 else
207 error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'curvature_tolerant'");
208 iarg_++;
209 hasargs = true;
210 } else if (strcmp(arg[iarg_], "extrude_planar") == 0) {
211 if (narg < iarg_+2)
212 error->fix_error(FLERR,this,"not enough arguments for 'extrude_planar'");
213 iarg_++;
214 extrude_mesh_ = true;
215 extrusion_length_ = force->numeric(FLERR,arg[iarg_]);
216 iarg_++;
217 hasargs = true;
218 } else if (strcmp(style,"mesh/surface") == 0) {
219 char *errmsg = new char[strlen(arg[iarg_])+50];
220 sprintf(errmsg,"unknown keyword or wrong keyword order: %s", arg[iarg_]);
221 error->fix_error(FLERR,this,errmsg);
222 delete []errmsg;
223 }
224 }
225
226 // get list of available mesh modules
227 std::map<std::string, MeshModuleCreator> *module_map = new std::map<std::string, MeshModuleCreator>();
228 std::map<std::string, std::vector<std::string> > *module_restrict_map = new std::map<std::string, std::vector<std::string> >();
229 #define MESHMODULE_CLASS
230 #define MeshModuleStyle(key, Class) \
231 (*module_map)[#key] = &meshmodule_creator<Class>;
232 #define MeshModuleRestrict(key1, key2) \
233 if (module_restrict_map->find(#key1) != module_restrict_map->end()) \
234 (*module_restrict_map)[#key1].push_back(#key2); \
235 else \
236 module_restrict_map->insert( std::pair<std::string, std::vector<std::string> >(#key1, std::vector<std::string>(1,#key2))); \
237 if (module_restrict_map->find(#key2) != module_restrict_map->end()) \
238 (*module_restrict_map)[#key2].push_back(#key1); \
239 else \
240 module_restrict_map->insert( std::pair<std::string, std::vector<std::string> >(#key2, std::vector<std::string>(1,#key1)));
241 #include "style_mesh_module.h"
242 #undef MeshModuleStyle
243 #undef MESHMODULE_CLASS
244
245 // parse style name of fix
246 std::string style(arg[2]);
247 std::istringstream data(style);
248 std::string line;
249 char separator = '/';
250 std::vector<std::string> module_list;
251
252 bool is_planar = false;
253
254 while (std::getline(data, line, separator))
255 {
256 if (line.compare("mesh") == 0 || line.compare("surface") == 0)
257 { /* do nothing */ }
258 else if (line.compare("planar") == 0) //support for the mesh/surface/planar style
259 {
260 if (active_mesh_modules.size() == 0)
261 is_planar = true;
262 else
263 error->fix_error(FLERR,this,"Fix mesh/surface contains the keyword planar and others. The planar keyword is only allowed on its own, i.e. mesh/surface/planar");
264 }
265 else if (module_map->find(line) != module_map->end()) // mesh modules
266 {
267 if (is_planar)
268 error->fix_error(FLERR,this,"Fix mesh/surface contains the keyword planar and others. The planar keyword is only allowed on its own, i.e. mesh/surface/planar");
269 MeshModuleCreator mm_creator = (*module_map)[line];
270 active_mesh_modules.insert(std::pair<std::string, MeshModule*>(line, mm_creator(lmp, iarg_, narg, arg, this)));
271 mesh_module_order.push_back(line);
272
273 }
274 else
275 {
276 char *errmsg = new char[strlen(arg[2])+50];
277 sprintf(errmsg,"unknown keyword in style definition: %s (%s)", arg[2], line.c_str());
278 error->fix_error(FLERR,this,errmsg);
279 delete []errmsg;
280 }
281 }
282
283 // check if we have illegal module combinations
284 std::map<std::string, MeshModule*>::iterator it;
285 for(it = active_mesh_modules.begin(); it != active_mesh_modules.end(); it++)
286 {
287 // check if the name of this module appears in the restrict map
288 std::map<std::string, std::vector<std::string> >::iterator cur_restrict_map = module_restrict_map->find(it->first);
289 if (cur_restrict_map != module_restrict_map->end())
290 {
291 // if yes check the vector
292 std::vector<std::string> *cur_restrict = &(cur_restrict_map->second);
293 std::map<std::string, MeshModule*>::iterator it2 = it;
294 // loop over remaining modlues and check if they are in the vector
295 for (std::advance(it2, 1); it2 != active_mesh_modules.end(); it2++)
296 {
297 std::vector<std::string>::iterator other_module = std::find(cur_restrict->begin(), cur_restrict->end(), it2->first);
298 if (other_module != cur_restrict->end()) {
299 char errormsg[200];
300 sprintf(errormsg, "Conflicting mesh modules loaded: %s and %s\n", cur_restrict_map->first.c_str(), other_module->c_str());
301 error->fix_error(FLERR, this, errormsg);
302 }
303 }
304 }
305 }
306
307 // check if these modules want to write vectors
308 for(it = active_mesh_modules.begin(); it != active_mesh_modules.end(); it++)
309 size_vector += it->second->get_num_vector_components();
310 if (size_vector > 0)
311 vector_flag = 1;
312
313 delete module_map;
314 delete module_restrict_map;
315
316 if (iarg_ < narg) {
317 char *errmsg = new char[strlen(arg[iarg_])+50];
318 sprintf(errmsg,"unknown keyword or wrong keyword order: %s", arg[iarg_]);
319 error->fix_error(FLERR,this,errmsg);
320 delete []errmsg;
321 }
322 }
323
324 /* ---------------------------------------------------------------------- */
325
~FixMeshSurface()326 FixMeshSurface::~FixMeshSurface()
327 {
328 delete [] vSurfStrX_;
329 delete [] vSurfStrY_;
330 delete [] vSurfStrZ_;
331 delete [] omegaStr_;
332 if (extrusion_tri_nodes_)
333 delete [] extrusion_tri_nodes_;
334 }
335
336 /* ---------------------------------------------------------------------- */
337
post_create()338 void FixMeshSurface::post_create()
339 {
340 FixMesh::post_create();
341
342 if(curvature_ > 0.)
343 triMesh()->setCurvature(curvature_);
344
345 if(curvature_tolerant_)
346 triMesh()->setCurvatureTolerant(curvature_tolerant_);
347
348 if(velFlag_ && angVelFlag_)
349 error->fix_error(FLERR,this,"cannot use 'surface_vel' and 'surface_ang_vel' together");
350
351 if(velFlag_)
352 initVel();
353
354 if(angVelFlag_)
355 initAngVel();
356
357 std::vector<std::string>::iterator it;
358 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
359 active_mesh_modules[*it]->post_create();
360 }
361
362 /* ---------------------------------------------------------------------- */
363
post_create_pre_restart()364 void FixMeshSurface::post_create_pre_restart()
365 {
366 std::vector<std::string>::iterator it;
367 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
368 active_mesh_modules[*it]->post_create_pre_restart();
369 }
370
371 /* ---------------------------------------------------------------------- */
372
init()373 void FixMeshSurface::init()
374 {
375 FixMesh::init();
376
377 if (vSurfStrX_)
378 {
379 vSurfVarX_ = input->variable->find(vSurfStrX_);
380 if (vSurfVarX_ < 0)
381 error->all(FLERR, "Variable name for fix mesh/surface surfaceVelX does not exist");
382 if (input->variable->equalstyle(vSurfVarX_))
383 vSurfStyleX_ = EQUAL;
384 else
385 error->all(FLERR, "Variable for fix mesh/surface surfaceVelX has invalid style");
386 }
387 if (vSurfStrY_)
388 {
389 vSurfVarY_ = input->variable->find(vSurfStrY_);
390 if (vSurfVarY_ < 0)
391 error->all(FLERR, "Variable name for fix mesh/surface surfaceVelY does not exist");
392 if (input->variable->equalstyle(vSurfVarY_))
393 vSurfStyleY_ = EQUAL;
394 else
395 error->all(FLERR, "Variable for fix mesh/surface surfaceVelY has invalid style");
396 }
397 if (vSurfStrZ_)
398 {
399 vSurfVarZ_ = input->variable->find(vSurfStrZ_);
400 if (vSurfVarZ_ < 0)
401 error->all(FLERR, "Variable name for fix mesh/surface surfaceVelZ does not exist");
402 if (input->variable->equalstyle(vSurfVarZ_))
403 vSurfStyleZ_ = EQUAL;
404 else
405 error->all(FLERR, "Variable for fix mesh/surface surfaceVelZ has invalid style");
406 }
407
408 if (omegaStr_)
409 {
410 omegaVar_ = input->variable->find(omegaStr_);
411 if (omegaVar_ < 0)
412 error->all(FLERR, "Variable name for fix mesh/surface omega does not exist");
413 if (input->variable->equalstyle(omegaVar_))
414 omegaStyle_ = EQUAL;
415 else
416 error->all(FLERR, "Variable for fix mesh/surface omega has invalid style");
417 }
418
419 std::vector<std::string>::iterator it;
420 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
421 active_mesh_modules[*it]->init();
422 }
423
424 /* ---------------------------------------------------------------------- */
425
setup(int vflag)426 void FixMeshSurface::setup(int vflag)
427 {
428 FixMesh::setup(vflag);
429
430 std::vector<std::string>::iterator it;
431 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
432 active_mesh_modules[*it]->setup(vflag);
433 }
434
435 /* ---------------------------------------------------------------------- */
436
pre_delete(bool unfixflag)437 void FixMeshSurface::pre_delete(bool unfixflag)
438 {
439 if(unfixflag && n_dump_active_ > 0)
440 error->fix_error(FLERR,this,"can not unfix while dump command is active on mesh");
441
442 if(unfixflag && fix_contact_history_mesh_)
443 error->fix_error(FLERR,this,"can not unfix while fix wall/gran command is active on mesh; need to unfix fix wall/gran first, then mesh");
444
445 FixMesh::pre_delete(unfixflag);
446
447 // contact tracker and neighlist are created via fix wall/gran
448 deleteWallNeighList();
449 deleteAllOtherNeighList();
450 deleteContactHistory();
451 }
452
453 /* ---------------------------------------------------------------------- */
454
setmask()455 int FixMeshSurface::setmask()
456 {
457 int mask = FixMesh::setmask();
458
459 if(velFlag_ && (vSurfVarX_ >= 0 || vSurfVarY_ >= 0 || vSurfVarZ_ >= 0))
460 mask |= PRE_FORCE;
461
462 if(angVelFlag_ && omegaVar_ >= 0)
463 mask |= PRE_FORCE;
464
465 std::vector<std::string>::iterator it;
466 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
467 mask |= active_mesh_modules[*it]->setmask();
468
469 return mask;
470 }
471
472 /* ---------------------------------------------------------------------- */
473
setup_pre_force(int vflag)474 void FixMeshSurface::setup_pre_force(int vflag)
475 {
476 FixMesh::setup_pre_force(vflag);
477
478 // create neigh list for owned and local elements
479
480 if(meshNeighlist())
481 meshNeighlist()->initializeNeighlist();
482
483 std::vector<std::string>::iterator it;
484 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
485 active_mesh_modules[*it]->setup_pre_force(vflag);
486
487 if (extrude_mesh_ && !extrusion_created_)
488 {
489 if (update->nsteps > 0)
490 error->all(FLERR, "Extrude mesh requires a run 0 after its definition");
491 mesh()->extrudePlanarMesh(extrusion_length_, extrusion_tri_nodes_, extrusion_tri_count_);
492 extrusion_created_ = true;
493 can_create_mesh_ = true;
494 }
495
496 }
497
498 /* ---------------------------------------------------------------------- */
499
pre_force(int vflag)500 void FixMeshSurface::pre_force(int vflag)
501 {
502 FixMesh::pre_force(vflag);
503
504 if(velFlag_ && (vSurfVarX_ >= 0 || vSurfVarY_ >= 0 || vSurfVarZ_ >= 0))
505 setVel();
506
507 if(angVelFlag_ && omegaVar_ >= 0)
508 setAngVel();
509
510 std::vector<std::string>::iterator it;
511 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
512 active_mesh_modules[*it]->pre_force(vflag);
513 }
514
515 /* ----------------------------------------------------------------------
516 initial integrate for mesh modules
517 ------------------------------------------------------------------------- */
518
initial_integrate(int vflag)519 void FixMeshSurface::initial_integrate(int vflag)
520 {
521 std::vector<std::string>::iterator it;
522 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
523 active_mesh_modules[*it]->initial_integrate(vflag);
524 }
525
526 /* ----------------------------------------------------------------------
527 reverse comm for mesh
528 ------------------------------------------------------------------------- */
529
final_integrate()530 void FixMeshSurface::final_integrate()
531 {
532 // Mesh modules final integrate before reverse communication of properties
533 std::vector<std::string>::iterator it;
534 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
535 active_mesh_modules[*it]->final_integrate_pre_comm();
536
537 // Reverse communication of mesh properties happens in here
538 FixMesh::final_integrate();
539
540 // Default final integrate of mesh modules after reverse communication
541 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
542 active_mesh_modules[*it]->final_integrate();
543 }
544
545 /* ----------------------------------------------------------------------
546 End of step call to modules
547 ------------------------------------------------------------------------- */
548
end_of_step()549 void FixMeshSurface::end_of_step()
550 {
551 std::vector<std::string>::iterator it;
552 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
553 active_mesh_modules[*it]->end_of_step();
554 }
555
556 /* ----------------------------------------------------------------------
557 Checks if any modules compute a vector
558 ------------------------------------------------------------------------- */
559
compute_vector(int n)560 double FixMeshSurface::compute_vector(int n)
561 {
562 std::vector<std::string>::iterator it;
563 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
564 {
565 const int num = active_mesh_modules[*it]->get_num_vector_components();
566 if (n < num)
567 return active_mesh_modules[*it]->compute_vector(n);
568 else
569 n -= num;
570 }
571 if (n == 0)
572 error->fix_error(FLERR, this, "Internal error");
573
574 return 0.0;
575 }
576
577 /* ----------------------------------------------------------------------
578 called from fix wall/gran out of post_create()
579 ------------------------------------------------------------------------- */
580
createWallNeighList(int igrp)581 void FixMeshSurface::createWallNeighList(int igrp)
582 {
583 if(fix_mesh_neighlist_) return;
584 char *neighlist_name = new char[strlen(id)+1+20];
585 sprintf(neighlist_name,"wall_neighlist_%s",id);
586
587 const char *fixarg[4];
588 fixarg[0]= neighlist_name;
589 fixarg[1]= "all";
590 fixarg[2]= "neighlist/mesh";
591 fixarg[3]= id;
592 modify->add_fix(4,const_cast<char**>(fixarg));
593
594 fix_mesh_neighlist_ =
595 static_cast<FixNeighlistMesh*>(modify->find_fix_id(neighlist_name));
596
597 // fix added with "all", correct this now
598 // groupbit identical to groupbit of this fix
599 fix_mesh_neighlist_->igroup = igroup;
600 fix_mesh_neighlist_->groupbit = group->bitmask[igroup];
601
602 // neighlist build will use this: merging of wall and mesh groupbit
603 fix_mesh_neighlist_->groupbit_wall_mesh = group->bitmask[igroup] | igrp;
604
605 delete []neighlist_name;
606
607 /*
608
609 fix_mesh_neighlist_->pre_neighbor();
610 fix_mesh_neighlist_->pre_force(0);
611 */
612 }
613
deleteWallNeighList()614 void FixMeshSurface::deleteWallNeighList()
615 {
616
617 if(fix_mesh_neighlist_) {
618 modify->delete_fix(fix_mesh_neighlist_->id,true);
619 fix_mesh_neighlist_ = NULL;
620 }
621 }
622
623 /* ----------------------------------------------------------------------
624 called from fix massflow/mesh out of post_create()
625 ------------------------------------------------------------------------- */
626
createOtherNeighList(int igrp,const char * nId)627 class FixNeighlistMesh* FixMeshSurface::createOtherNeighList(int igrp,const char *nId)
628 {
629 FixNeighlistMesh* neighlist;
630
631 char *neighlist_name = new char[strlen(id)+1+20+strlen(nId)+1];
632 sprintf(neighlist_name,"neighlist_%s_%s",nId,id);
633
634 if(modify->find_fix_id(neighlist_name))
635 error->fix_error(FLERR,this,"must not use the same mesh for fix massflow/mesh with same group");
636
637 const char *fixarg[5];
638 fixarg[0]= neighlist_name;
639 fixarg[1]= "all";
640 fixarg[2]= "neighlist/mesh";
641 fixarg[3]= id;
642 fixarg[4]= "other_yes";
643 modify->add_fix(5,const_cast<char**>(fixarg));
644
645 neighlist =
646 static_cast<FixNeighlistMesh*>(modify->find_fix_id(neighlist_name));
647
648 // fix added with "all", correct this now
649 neighlist->igroup = igrp;
650 neighlist->groupbit = group->bitmask[igrp];
651
652 list_other_neighlist_[neighlist_name] = neighlist;
653
654 delete []neighlist_name;
655
656 return neighlist;
657 }
658
deleteOtherNeighList(const char * nId)659 void FixMeshSurface::deleteOtherNeighList(const char *nId)
660 {
661 char *neighlist_name = new char[strlen(id)+1+20+strlen(nId)+1];
662 sprintf(neighlist_name,"neighlist_%s_%s",nId,id);
663
664 std::map<std::string,FixNeighlistMesh*>::iterator it = list_other_neighlist_.find(neighlist_name);
665 if (it != list_other_neighlist_.end()) {
666 modify->delete_fix(it->second->id,true);
667 list_other_neighlist_.erase(it);
668 }
669
670 delete []neighlist_name;
671 }
672
deleteAllOtherNeighList()673 void FixMeshSurface::deleteAllOtherNeighList()
674 {
675 for ( std::map<std::string,FixNeighlistMesh*>::iterator it = list_other_neighlist_.begin();
676 it != list_other_neighlist_.end();
677 ++ it) {
678 modify->delete_fix(it->second->id,true);
679 }
680 list_other_neighlist_.clear();
681 }
682
683 /* ----------------------------------------------------------------------
684 called from fix wall/gran out of post_create()
685 ------------------------------------------------------------------------- */
686
createContactHistory(int dnum)687 void FixMeshSurface::createContactHistory(int dnum)
688 {
689 if(fix_contact_history_mesh_) return;
690
691 // create a contact tracker for the mesh
692 char *contacthist_name = new char[strlen(id)+1+8];
693 char *my_id = new char[strlen(id)+1];
694 sprintf(contacthist_name,"tracker_%s",id);
695 sprintf(my_id,"%s",id);
696
697 char dnum_char[10];
698 sprintf(dnum_char,"%d",dnum);
699
700 const char *fixarg[5];
701 fixarg[0] = contacthist_name;
702 fixarg[1] = "all";
703 fixarg[2] = "contacthistory/mesh";
704 fixarg[3] = dnum_char;
705 fixarg[4] = my_id;
706
707 modify->add_fix(5,const_cast<char**>(fixarg));
708
709 fix_contact_history_mesh_ = static_cast<FixContactHistoryMesh*>(modify->find_fix_id(contacthist_name));
710
711 delete []contacthist_name;
712 delete []my_id;
713 }
714
deleteContactHistory()715 void FixMeshSurface::deleteContactHistory()
716 {
717 // contact tracker and neighlist are created via fix wall/gran
718
719 if(fix_contact_history_mesh_) {
720 modify->delete_fix(fix_contact_history_mesh_->id,true);
721 fix_contact_history_mesh_ = NULL;
722 }
723 }
724
725 /* ----------------------------------------------------------------------
726 called from fix wall/gran out of post_create()
727 ------------------------------------------------------------------------- */
728
createMeshforceContact()729 void FixMeshSurface::createMeshforceContact()
730 {
731 if(fix_meshforce_contact_) return;
732
733 const char *fixarg[19];
734 char fixid[200],propertyid[200];
735 sprintf(fixid,"contactforces_%s",id);
736 sprintf(propertyid,"contactforces_%s",id);
737 fixarg[0]=fixid;
738 fixarg[1]="all";
739 fixarg[2]="contactproperty/atom/wall";
740 fixarg[3]=propertyid;
741 fixarg[4]="6";
742 fixarg[5]="fx";
743 fixarg[6]="0";
744 fixarg[7]="fy";
745 fixarg[8]="0";
746 fixarg[9]="fz";
747 fixarg[10]="0";
748 fixarg[11]="tx";
749 fixarg[12]="0";
750 fixarg[13]="ty";
751 fixarg[14]="0";
752 fixarg[15]="tz";
753 fixarg[16]="0";
754 fixarg[17]="mesh";
755 fixarg[18]=this->id;
756 modify->add_fix(19,const_cast<char**>(fixarg));
757 fix_meshforce_contact_ = static_cast<FixContactPropertyAtomWall*>(modify->find_fix_id(fixid));
758 }
759
createMeshforceContactStress()760 void FixMeshSurface::createMeshforceContactStress()
761 {
762 if(fix_meshforce_contact_stress_) return;
763
764 const char *fixarg[25];
765 char fixid[200],propertyid[200];
766 sprintf(fixid,"contactforces_stress_%s",id);
767 sprintf(propertyid,"contactforces_stress_%s",id);
768 fixarg[0]=fixid;
769 fixarg[1]="all";
770 fixarg[2]="contactproperty/atom/wall";
771 fixarg[3]=propertyid;
772 fixarg[4]="9";
773 fixarg[5]="fx";
774 fixarg[6]="0";
775 fixarg[7]="fy";
776 fixarg[8]="0";
777 fixarg[9]="fz";
778 fixarg[10]="0";
779 fixarg[11]="deltax";
780 fixarg[12]="0";
781 fixarg[13]="deltay";
782 fixarg[14]="0";
783 fixarg[15]="deltaz";
784 fixarg[16]="0";
785 fixarg[17]="vx";
786 fixarg[18]="0";
787 fixarg[19]="vy";
788 fixarg[20]="0";
789 fixarg[21]="vz";
790 fixarg[22]="0";
791 fixarg[23]="mesh";
792 fixarg[24]=this->id;
793 modify->add_fix(25,const_cast<char**>(fixarg));
794 fix_meshforce_contact_stress_ = static_cast<FixContactPropertyAtomWall*>(modify->find_fix_id(fixid));
795 }
796
createMulticontactData()797 void FixMeshSurface::createMulticontactData()
798 {
799 if(fix_mesh_multicontact_data_) return;
800
801 // create a new per contact property which will contain the data for the computation according to Brodu et. al. 2016
802 // surfPosIJ will contain the position of the contact surface ij, realtive to position i
803 // normalForce will contain the normal component of the contact force
804 const char *fixarg[17];
805 char fixid[200];
806 sprintf(fixid,"multicontactData_%s",id);
807 fixarg[0]=fixid;
808 fixarg[1]="all";
809 fixarg[2]="contactproperty/atom/wall";
810 fixarg[3]=fixid;
811 fixarg[4]="4";
812 fixarg[5]="surfPosIJ_x";
813 fixarg[6]="0";
814 fixarg[7]="surfPosIJ_y";
815 fixarg[8]="0";
816 fixarg[9]="surfPosIJ_z";
817 fixarg[10]="0";
818 fixarg[11]="normalForce";
819 fixarg[12]="0";
820 fixarg[13]="mesh";
821 fixarg[14]=id;
822 fixarg[15]="reset";
823 fixarg[16]="no";
824 modify->add_fix(17,const_cast<char**>(fixarg));
825 fix_mesh_multicontact_data_ = static_cast<FixContactPropertyAtomWall*>(modify->find_fix_id(fixid));
826 }
827
deleteMeshforceContact()828 void FixMeshSurface::deleteMeshforceContact()
829 {
830 // contact tracker and neighlist are created via fix wall/gran
831
832 if(fix_meshforce_contact_) {
833 modify->delete_fix(fix_meshforce_contact_->id,true);
834 fix_meshforce_contact_ = NULL;
835 }
836 }
837
deleteMeshforceContactStress()838 void FixMeshSurface::deleteMeshforceContactStress()
839 {
840 // contact tracker and neighlist are created via fix wall/gran
841
842 if(fix_meshforce_contact_stress_) {
843 modify->delete_fix(fix_meshforce_contact_stress_->id,true);
844 fix_meshforce_contact_stress_ = NULL;
845 }
846 }
847
deleteMeshMulticontactData()848 void FixMeshSurface::deleteMeshMulticontactData()
849 {
850 // contact tracker and neighlist are created via fix wall/gran
851
852 if(fix_mesh_multicontact_data_) {
853 modify->delete_fix(fix_mesh_multicontact_data_->id,true);
854 fix_mesh_multicontact_data_ = NULL;
855 }
856 }
857
858 /* ----------------------------------------------------------------------
859 sets mesh velocity for conveyor model
860 ------------------------------------------------------------------------- */
861
initVel()862 void FixMeshSurface::initVel()
863 {
864 if (vSurfStyleX_ == EQUAL)
865 {
866 modify->clearstep_compute();
867 vSurf_[0] = input->variable->compute_equal(vSurfVarX_);
868 modify->addstep_compute(update->ntimestep + 1);
869 }
870 if (vSurfStyleY_ == EQUAL)
871 {
872 modify->clearstep_compute();
873 vSurf_[1] = input->variable->compute_equal(vSurfVarY_);
874 modify->addstep_compute(update->ntimestep + 1);
875 }
876 if (vSurfStyleZ_ == EQUAL)
877 {
878 modify->clearstep_compute();
879 vSurf_[2] = input->variable->compute_equal(vSurfVarZ_);
880 modify->addstep_compute(update->ntimestep + 1);
881 }
882 double conv_vel[3];
883 vectorCopy3D(vSurf_,conv_vel);
884
885 // register new mesh property v
886 mesh()->prop().addGlobalProperty< VectorContainer<double,3> >("v","comm_exchange_borders","frame_invariant","restart_no");
887 mesh()->prop().setGlobalProperty< VectorContainer<double,3> >("v",conv_vel);
888
889 // register new element property v - error if exists already via moving mesh
890 mesh()->prop().addElementProperty<MultiVectorContainer<double,3,3> >("v","comm_exchange_borders","frame_invariant","restart_no");
891
892 setVel();
893 }
894
setVel()895 void FixMeshSurface::setVel()
896 {
897 if (vSurfStyleX_ == EQUAL)
898 {
899 modify->clearstep_compute();
900 vSurf_[0] = input->variable->compute_equal(vSurfVarX_);
901 modify->addstep_compute(update->ntimestep + 1);
902 }
903 if (vSurfStyleY_ == EQUAL)
904 {
905 modify->clearstep_compute();
906 vSurf_[1] = input->variable->compute_equal(vSurfVarY_);
907 modify->addstep_compute(update->ntimestep + 1);
908 }
909 if (vSurfStyleZ_ == EQUAL)
910 {
911 modify->clearstep_compute();
912 vSurf_[2] = input->variable->compute_equal(vSurfVarZ_);
913 modify->addstep_compute(update->ntimestep + 1);
914 }
915 double conv_vel[3];
916 int size, nVec;
917 double scp, tmp[3], facenormal[3], ***v_node;
918 vectorCopy3D(vSurf_,conv_vel);
919 double conv_vSurf_mag = vectorMag3D(conv_vel);
920
921 size = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->size();
922 nVec = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->nVec();
923
924 v_node = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->begin();
925
926 // set mesh velocity
927 TriMesh *trimesh = triMesh();
928
929 for (int i = 0; i < size; i++)
930 {
931 trimesh->surfaceNorm(i,facenormal);
932 scp = vectorDot3D(conv_vel,facenormal);
933 vectorScalarMult3D(facenormal,scp,tmp);
934 for(int j = 0; j < nVec; j++)
935 {
936 vectorSubtract3D(conv_vel,tmp,v_node[i][j]);
937 if(vectorMag3D(v_node[i][j]) > 0.)
938 {
939 vectorScalarDiv3D(v_node[i][j],vectorMag3D(v_node[i][j]));
940 vectorScalarMult3D(v_node[i][j],conv_vSurf_mag);
941
942 }
943 }
944 }
945 }
946
947 /* ----------------------------------------------------------------------
948 sets mesh velocity for rotation model
949 ------------------------------------------------------------------------- */
950
initAngVel()951 void FixMeshSurface::initAngVel()
952 {
953 // register new element property v - error if exists already via moving mesh
954 mesh()->prop().addElementProperty<MultiVectorContainer<double,3,3> >("v","comm_exchange_borders","frame_invariant","restart_no");
955
956 setAngVel();
957 }
958
setAngVel()959 void FixMeshSurface::setAngVel()
960 {
961 if (omegaStyle_ == EQUAL)
962 {
963 modify->clearstep_compute();
964 omegaSurf_ = input->variable->compute_equal(omegaVar_);
965 modify->addstep_compute(update->ntimestep + 1);
966 }
967 int size, nVec;
968 double rot_origin[3],rot_axis[3],rot_omega;
969 vectorCopy3D(origin_,rot_origin);
970 vectorCopy3D(axis_,rot_axis);
971 rot_omega = omegaSurf_;
972 double tmp[3], scp, unitAxis[3], tangComp[3], Utang[3], surfaceV[3];
973 double node[3],facenormal[3], magAxis, magUtang, ***v_node;
974
975 size = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->size();
976 nVec = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->nVec();
977
978 v_node = mesh()->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->begin();
979
980 // calculate unit vector of rotation axis
981 magAxis = vectorMag3D(rot_axis);
982 vectorScalarDiv3D(rot_axis,magAxis,unitAxis);
983
984 TriMesh *trimesh = triMesh();
985 for (int i = 0; i < size; i++)
986 {
987 trimesh->surfaceNorm(i,facenormal);
988 // number of nodes per face (3)
989 for(int j = 0; j < nVec; j++)
990 {
991 // position of node - origin of rotation (to get lever arm)
992 trimesh->node(i,j,node);
993 vectorSubtract3D(node,rot_origin,surfaceV);
994
995 // lever arm X rotational axis = tangential component
996 vectorCross3D(surfaceV,unitAxis,tangComp);
997
998 // multiplying by omega scales the tangential component to give tangential velocity
999 vectorScalarMult3D(tangComp,-rot_omega,Utang);
1000
1001 // EPSILON is wall velocity, not rotational omega
1002 //if(vectorMag3D(Utang) < EPSILON_V)
1003 // error->fix_error(FLERR,this,"Rotation velocity too low");
1004
1005 scp = vectorDot3D(Utang, facenormal);
1006 vectorScalarMult3D(facenormal,scp,tmp);
1007
1008 // removes components normal to wall
1009 vectorSubtract3D(Utang,tmp,v_node[i][j]);
1010 magUtang = vectorMag3D(Utang);
1011
1012 if(vectorMag3D(v_node[i][j]) > 0.)
1013 {
1014 vectorScalarDiv3D(v_node[i][j],vectorMag3D(v_node[i][j]));
1015 vectorScalarMult3D(v_node[i][j],magUtang);
1016 }
1017 }
1018 }
1019 }
1020
1021 /* ----------------------------------------------------------------------
1022 Returns a pointer to a specific mesh module if it is loaded
1023 ------------------------------------------------------------------------- */
1024
get_module(std::string name)1025 MeshModule* FixMeshSurface::get_module(std::string name)
1026 {
1027 std::map<std::string, MeshModule*>::iterator it = active_mesh_modules.find(name);
1028 if (it != active_mesh_modules.end())
1029 return it->second;
1030 else
1031 return NULL;
1032 }
1033
1034 /* ----------------------------------------------------------------------
1035 Creates a new mesh module
1036 ------------------------------------------------------------------------- */
1037
1038 template <typename T>
meshmodule_creator(LAMMPS * lmp,int & iarg_,int narg,char ** arg,FixMeshSurface * fix_mesh)1039 MeshModule *FixMeshSurface::meshmodule_creator(LAMMPS *lmp, int &iarg_, int narg, char **arg, FixMeshSurface *fix_mesh)
1040 {
1041 return new T(lmp, iarg_, narg, arg, fix_mesh);
1042 }
1043
add_particle_contribution(int ip,double * frc,double * delta,int iTri,double * v_wall)1044 void FixMeshSurface::add_particle_contribution(int ip, double *frc, double *delta, int iTri, double *v_wall)
1045 {
1046 std::vector<std::string>::iterator it;
1047 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
1048 active_mesh_modules[*it]->add_particle_contribution(ip, frc, delta, iTri, v_wall);
1049 }
1050
1051 /* ----------------------------------------------------------------------
1052 Returns if the mesh tracks the stress.
1053 This function should not used within a performance critical code.
1054 Instead, save the value during an init/setup function.
1055 ------------------------------------------------------------------------- */
1056
trackStress()1057 bool FixMeshSurface::trackStress()
1058 {
1059 MeshModuleStress *mms = dynamic_cast<MeshModuleStress*>(get_module("stress"));
1060 if (mms)
1061 return mms->trackStress();
1062 else
1063 return false;
1064 }
1065
modify_param(int narg,char ** arg)1066 int FixMeshSurface::modify_param(int narg, char **arg)
1067 {
1068 int ret = 0;
1069 std::string module(arg[0]);
1070 std::size_t slash = module.find_first_of('/');
1071
1072 if (slash != std::string::npos)
1073 {
1074 MeshModule *mm_modify = get_module(module.substr(0,slash));
1075 if (!mm_modify)
1076 {
1077 std::string errmsg = "Could not find appropriate mesh module \"";
1078 errmsg.append(module.substr(0,slash));
1079 errmsg.append("\" in modify_param");
1080 error->fix_error(FLERR, this, errmsg.c_str());
1081 }
1082 ret = mm_modify->modify_param(narg, arg);
1083 return ret;
1084 }
1085 else
1086 {
1087 // Basic divide and conquer algorithm for modify_param in mesh modules. Note that this can cause issues if
1088 // two mesh modules use the same keyword.
1089 std::vector<std::string>::iterator it;
1090 for(it = mesh_module_order.begin(); it != mesh_module_order.end(); it++)
1091 {
1092 ret = active_mesh_modules[*it]->modify_param(narg, arg);
1093 if (ret > 0)
1094 {
1095 std::string warnmsg = "Using deprecated modify_param for attribute \"";
1096 warnmsg.append(module);
1097 warnmsg.append("\". Consider using ");
1098 warnmsg.append(*it);
1099 warnmsg.append("/");
1100 warnmsg.append(module);
1101 warnmsg.append(" in your input file to explicitly specify the mesh module the argument is used for.\n");
1102 error->warning(FLERR, warnmsg.c_str());
1103 return ret;
1104 }
1105 }
1106 }
1107 return ret;
1108 }
1109