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