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 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #include <cmath>
43 #include <stdlib.h>
44 #include <string.h>
45 #include "atom.h"
46 #include "atom_vec.h"
47 #include "comm.h"
48 #include "modify.h"
49 #include "force.h"
50 #include "memory.h"
51 #include "update.h"
52 #include "error.h"
53 #include "fix_mesh_surface.h"
54 #include "fix_neighlist_mesh.h"
55 #include "fix_multisphere.h"
56 #include "fix_property_atom.h"
57 #include "mpi_liggghts.h"
58 #include "math_extra_liggghts.h"
59 #include "fix_massflow_mesh.h"
60
61 #include <algorithm>
62
63 using namespace LAMMPS_NS;
64 using namespace MathExtraLiggghts;
65 using namespace FixConst;
66
67 /* ---------------------------------------------------------------------- */
68
FixMassflowMesh(LAMMPS * lmp,int narg,char ** arg)69 FixMassflowMesh::FixMassflowMesh(LAMMPS *lmp, int narg, char **arg) :
70 Fix(lmp, narg, arg),
71 delete_atoms_(false),
72 mass_deleted_(0.),
73 nparticles_deleted_(0),
74 once_(false),
75 fix_orientation_(0),
76 fix_counter_(0),
77 fix_mesh_(0),
78 fix_neighlist_(0),
79 fix_volumeweight_ms_(0),
80 havePointAtOutlet_(false),
81 insideOut_(false),
82 mass_(0.),
83 nparticles_(0.),
84 fix_property_(0),
85 property_sum_(0.),
86 screenflag_(false),
87 fp_(0),
88 writeTime_(false),
89 mass_last_(0.),
90 nparticles_last_(0.),
91 t_count_(0.),
92 delta_t_(0.),
93 reset_t_count_(true)
94 {
95 vectorZeroize3D(nvec_);
96 vectorZeroize3D(pref_);
97 vectorZeroize3D(sidevec_);
98 vectorZeroize3D(pointAtOutlet_);
99
100 // parse args for this class
101
102 iarg_ = 3;
103
104 bool hasargs = true;
105 while(iarg_ < narg && hasargs)
106 {
107 hasargs = false;
108
109 if(strcmp(arg[iarg_],"vec_side") == 0) {
110 if(narg < iarg_+4)
111 error->fix_error(FLERR,this,"not enough arguments for 'vec_side'");
112 iarg_++;
113 sidevec_[0] = atof(arg[iarg_++]);
114 sidevec_[1] = atof(arg[iarg_++]);
115 sidevec_[2] = atof(arg[iarg_++]);
116 if(vectorMag3D(sidevec_) == 0.)
117 error->fix_error(FLERR,this,"vec_side > 0 required");
118 hasargs = true;
119 } else if(strcmp(arg[iarg_],"mesh") == 0) {
120 if(narg < iarg_+2)
121 error->fix_error(FLERR,this,"not enough arguments for 'insert_stream'");
122 iarg_++;
123 fix_mesh_ = static_cast<FixMeshSurface*>(modify->find_fix_id_style(arg[iarg_++],"mesh/surface"));
124 if(!fix_mesh_)
125 error->fix_error(FLERR,this,"fix mesh ID does not exist");
126 hasargs = true;
127 } else if(strcmp(arg[iarg_],"sum_property") == 0) {
128 if(narg < iarg_+2)
129 error->fix_error(FLERR,this,"not enough arguments for 'sum_property'");
130 iarg_++;
131 fix_property_ = static_cast<FixPropertyAtom*>(modify->find_fix_property(arg[iarg_++],"property/atom","scalar",0,0,style));
132 hasargs = true;
133 } else if(strcmp(arg[iarg_],"count") == 0) {
134 if(narg < iarg_+2)
135 error->fix_error(FLERR,this,"not enough arguments for 'insert_stream'");
136 iarg_++;
137 if(strcmp(arg[iarg_],"once") == 0)
138 once_ = true;
139 else if(strcmp(arg[iarg_],"multiple") == 0)
140 once_ = false;
141 else
142 error->fix_error(FLERR,this,"expecting 'once' or 'multiple' after 'count'");
143 iarg_++;
144 hasargs = true;
145 } else if( strcmp(arg[iarg_],"writeTime") == 0) {
146 writeTime_ = true;
147 iarg_++;
148 hasargs = true;
149 } else if(strcmp(arg[iarg_],"point_at_outlet") == 0) {
150 if(narg < iarg_+4)
151 error->fix_error(FLERR,this,"not enough arguments for 'point_at_outlet'");
152 havePointAtOutlet_ = true;
153 iarg_++;
154 pointAtOutlet_[0] = atof(arg[iarg_++]);
155 pointAtOutlet_[1] = atof(arg[iarg_++]);
156 pointAtOutlet_[2] = atof(arg[iarg_++]);
157 hasargs = true;
158 } else if(strcmp(arg[iarg_],"inside_out") == 0) {
159 insideOut_ = true;
160 iarg_++;
161 if(!havePointAtOutlet_)
162 error->fix_error(FLERR,this,"the setting 'inside_out' has no meaning in case you do not use 'point_at_outlet'");
163 hasargs = true;
164 } else if (strcmp(arg[iarg_],"file") == 0 || strcmp(arg[iarg_],"append") == 0)
165 {
166 if(narg < iarg_+2)
167 error->fix_error(FLERR,this,"Illegal keyword entry");
168
169 char* filecurrent = new char[strlen(arg[iarg_+1]) + 8];
170 if (1 < comm->nprocs) //open a separate file for each processor
171 sprintf(filecurrent,"%s%s%d",arg[iarg_+1],".",comm->me);
172 else //open one file for proc 0
173 sprintf(filecurrent,"%s",arg[iarg_+1]);
174
175 if (strcmp(arg[iarg_],"file") == 0)
176 fp_ = fopen(filecurrent,"w");
177 else
178 fp_ = fopen(filecurrent,"a");
179 delete[] filecurrent;
180 if (fp_ == NULL) {
181 char str[512];
182 sprintf(str,"Cannot open file %s",arg[iarg_+1]);
183 error->fix_error(FLERR,this,str);
184 }
185 iarg_ += 2;
186 hasargs = true;
187 } else if (strcmp(arg[iarg_],"screen") == 0) {
188 if(narg < iarg_+2)
189 error->fix_error(FLERR,this,"Illegal keyword entry");
190 if (strcmp(arg[iarg_+1],"yes") == 0) screenflag_ = true;
191 else if (strcmp(arg[iarg_+1],"no") == 0) screenflag_ = false;
192 else error->all(FLERR,"Illegal fix print command");
193 iarg_ += 2;
194 hasargs = true;
195 } else if (strcmp(arg[iarg_],"delete_atoms") == 0) {
196 if(narg < iarg_+2)
197 error->fix_error(FLERR,this,"Illegal keyword entry");
198 if (strcmp(arg[iarg_+1],"yes") == 0) delete_atoms_ = true;
199 else if (strcmp(arg[iarg_+1],"no") == 0) delete_atoms_ = false;
200 else error->all(FLERR,"Illegal delete command");
201 iarg_ += 2;
202 hasargs = true;
203 } else if(strcmp(style,"massflow/mesh") == 0)
204 error->fix_error(FLERR,this,"unknown keyword");
205 }
206
207 if(fp_ && 1 < comm->nprocs && 0 == comm->me)
208 fprintf(screen,"**FixMassflowMesh: > 1 process - "
209 " will write to multiple files\n");
210 if(fp_)
211 {
212 //write header
213 fprintf(fp_,"# ID");
214
215 if(writeTime_)
216 fprintf(fp_," time ");
217
218 fprintf(fp_," diameter x y z u v w");
219
220 fprintf(fp_," (ex ey ez, color)\n");
221
222 fflush(fp_);
223 }
224
225 // error checks on necessary args
226
227 if(!once_ && delete_atoms_)
228 error->fix_error(FLERR,this,"using 'delete_atoms yes' requires 'count once'");
229 if( !once_ && havePointAtOutlet_)
230 error->fix_error(FLERR,this,"setting 'point_at_outlet' requires 'count once'");
231 if( vectorMag3D(sidevec_)==0. && !havePointAtOutlet_)
232 error->fix_error(FLERR,this,"expecting keyword 'vec_side'");
233 if (!fix_mesh_)
234 error->fix_error(FLERR,this,"expecting keyword 'mesh'");
235
236 // get reference point on face
237 // calculate normalvec
238 setRefPoint();
239
240 restart_global = 1;
241
242 vector_flag = 1;
243 size_vector = 6;
244 if(fix_property_)
245 size_vector = 7;
246 global_freq = 1;
247 }
248
249 /* ---------------------------------------------------------------------- */
250
~FixMassflowMesh()251 FixMassflowMesh::~FixMassflowMesh()
252 {
253 if(fp_)
254 fclose(fp_);
255 // do not delete fix_neighlist_, this will be handled by fix mesh/surface itself
256 fix_neighlist_ = NULL;
257 }
258
259 /* ---------------------------------------------------------------------- */
260
post_create()261 void FixMassflowMesh::post_create()
262 {
263 // add per-particle count flag
264
265 const char * fixarg[9];
266
267 sprintf(fixid_,"massflow_%s",id);
268 fixarg[0]=fixid_;
269 fixarg[1]="all";
270 fixarg[2]="property/atom";
271 fixarg[3]=fixid_;
272 fixarg[4]="scalar";
273 fixarg[5]="yes";
274 fixarg[6]="no";
275 fixarg[7]="no";
276 fixarg[8]="-1.";
277 modify->add_fix(9,const_cast<char**>(fixarg));
278
279 fix_counter_ = static_cast<FixPropertyAtom*>(modify->find_fix_property(fixid_,"property/atom","scalar",0,0,style));
280
281 // add neighbor list
282
283 fix_neighlist_ = fix_mesh_->createOtherNeighList(igroup,id);
284
285 fix_volumeweight_ms_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("volumeweight_ms","property/atom","scalar",0,0,style,false));
286
287 // need to find multisphere here to be able to add property
288
289 FixMultisphere *fix_ms = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
290 if(fix_ms)
291 {
292 Multisphere *ms = &fix_ms->data();
293 char property_name[200];
294 sprintf(property_name,"counter_ms_%s",id);
295
296 if(!ms->prop().getElementProperty< ScalarContainer<int> >(static_cast<const char*>(property_name)))
297 {
298 (ms->prop().addElementProperty< ScalarContainer<int> >(static_cast<const char*>(property_name),"comm_exchange_borders","frame_invariant", "restart_yes"))->setDefaultValue(-1);
299
300 }
301
302 if(delete_atoms_)
303 error->fix_error(FLERR,this,"can not use 'delete_atoms' with fix multisphere/*");
304 if(!once_)
305 error->fix_error(FLERR,this,"must use 'count once' with fix multisphere/*");
306 }
307
308 }
309
310 /* ---------------------------------------------------------------------- */
311
pre_delete(bool unfixflag)312 void FixMassflowMesh::pre_delete(bool unfixflag)
313 {
314 if (unfixflag) modify->delete_fix(fixid_);
315 }
316
317 /* ----------------------------------------------------------------------
318 initialize this fix
319 ------------------------------------------------------------------------- */
320
init()321 void FixMassflowMesh::init()
322 {
323
324 fix_volumeweight_ms_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("volumeweight_ms","property/atom","scalar",0,0,style,false));
325
326 if (atom->rmass_flag == 0)
327 error->fix_error(FLERR,this,"requires atoms have mass");
328
329 if(delete_atoms_ && 1 != atom->map_style)
330 error->fix_error(FLERR,this,"requires an atom map of type 'array', via an 'atom_modify map array' command");
331 /*
332 if(!fix_ms_ && static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0)))
333 error->fix_error(FLERR,this,"fix multisphere must come before fix massflow/mesh in input script");*/
334 }
335
336 /* ---------------------------------------------------------------------- */
337
setup(int vflag)338 void FixMassflowMesh::setup(int vflag)
339 {
340 // check if face planar
341
342 if(!fix_mesh_->triMesh()->isPlanar() && !havePointAtOutlet_)
343 error->fix_error(FLERR,this,"requires a planar face mass flow measurement or using 'point_at_outlet'");
344 }
345
346 /* ---------------------------------------------------------------------- */
347
setmask()348 int FixMassflowMesh::setmask()
349 {
350 int mask = 0;
351 mask |= POST_INTEGRATE;
352 if(delete_atoms_) mask |= PRE_EXCHANGE;
353 return mask;
354 }
355
356 /* ----------------------------------------------------------------------
357 evaluate mass which went thru face
358 all nearby particles
359 -1 = default value for counter, set if particle is not
360 a neighbor of the TriMesh
361 0 = was not on nvec_ last step and is thus eligible
362 1 = was on nvec_ side last step
363 2 = do not re-count, was counted already
364 ------------------------------------------------------------------------- */
365
post_integrate()366 void FixMassflowMesh::post_integrate()
367 {
368
369 int nlocal = atom->nlocal;
370 double **x = atom->x;
371 double **v = atom->v;
372 double *radius = atom->radius;
373 double *rmass = atom->rmass;
374 int *mask = atom->mask;
375 double *counter = fix_counter_->vector_atom;
376 double dot,delta[3]={}, bary[3];
377 double mass_this = 0.;
378 double nparticles_this = 0.;
379 double property_this = 0.;
380 double deltan;
381 int *tag = atom->tag;
382
383 class FixPropertyAtom* fix_color=static_cast<FixPropertyAtom*>(modify->find_fix_property("color","property/atom","scalar",0,0,style,false));
384 bool fixColFound = false;
385 if (fix_color) fixColFound=true;
386
387 fix_orientation_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("ex",
388 "property/atom","vector",0,0,style,false));
389
390 TriMesh *mesh = fix_mesh_->triMesh();
391 int nTriAll = mesh->sizeLocal() + mesh->sizeGhost();
392
393 // update reference point
394 if (fix_mesh_->triMesh()->isMoving() || fix_mesh_->triMesh()->isDeforming()) {
395 setRefPoint();
396 }
397
398 // update time for counter
399 // also store values for last invokation
400 t_count_ += update->dt;
401 if(!reset_t_count_)
402 {
403 nparticles_last_ = nparticles_;
404 mass_last_ = mass_;
405 reset_t_count_ = true;
406 }
407
408 // loop owned and ghost triangles
409 // count only if owned particle
410
411 for(int iTri = 0; iTri < nTriAll; iTri++)
412 {
413
414 const std::vector<int> & neighborList = fix_neighlist_->get_contact_list(iTri);
415 const int numneigh = neighborList.size();
416
417 for(int iNeigh = 0; iNeigh < numneigh; iNeigh++)
418 {
419 const int iPart = neighborList[iNeigh];
420
421 // skip ghost particles
422 if(iPart >= nlocal)
423 continue;
424
425 // skip particles not in fix group
426 if (!(mask[iPart] & groupbit))
427 continue;
428
429 // in case of once_ == true, ignore everything which has been already counted
430 if(compDouble(counter[iPart],2.))
431 continue;
432
433 int barySign;
434 deltan = fix_mesh_->triMesh()->resolveTriSphereContactBary(iPart,iTri,radius[iPart],x[iPart],delta,bary,barySign);
435
436 if(havePointAtOutlet_)
437 {
438 if(deltan < radius[iPart])
439 {
440 vectorSubtract3D(x[iPart],pointAtOutlet_,nvec_); //vector pointing to the particle location
441 dot = vectorDot3D(delta,nvec_);
442 }
443 else //particle is not overlapping with mesh, so continue
444 continue;
445
446 if(insideOut_) dot =-dot;
447 }
448 else
449 {
450 vectorSubtract3D(x[iPart],pref_,delta);
451 dot = vectorDot3D(delta,nvec_);
452 }
453
454 // first-time, just set 0 or 1 depending on what side of the mesh
455 if(compDouble(counter[iPart],-1.))
456 {
457 counter[iPart] = (dot <= 0.) ? 0. : 1.;
458
459 continue;
460 }
461
462 // particle is now on nvec_ side
463 if(dot > 0. && 7 == barySign)
464 {
465 //particle was not on nvec_ side before
466 if((compDouble(counter[iPart],0.)) ) // compDouble(counter[iPart],0.))
467 {
468
469 mass_this += (fix_volumeweight_ms_ ? fix_volumeweight_ms_->vector_atom[iPart] : 1.) * rmass[iPart];
470 nparticles_this += (fix_volumeweight_ms_ ? fix_volumeweight_ms_->vector_atom[iPart] : 1.);
471
472 if(fix_property_)
473 {
474
475 property_this += fix_property_->vector_atom[iPart];
476 }
477
478 if(delete_atoms_)
479 {
480 //reset counter to avoid problems with other fixes & mark to be deleted
481 counter[iPart] = -1.0;
482 atom_tags_delete_.push_back(atom->tag[iPart]);
483 }
484
485 if (screenflag_ && screen)
486 fprintf(screen," %d %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g \n ",
487 tag[iPart],2.*radius[iPart]/force->cg(atom->type[iPart]),
488 x[iPart][0],x[iPart][1],x[iPart][2],
489 v[iPart][0],v[iPart][1],v[iPart][2]);
490 if(fp_)
491 {
492 fprintf(fp_,"%d", tag[iPart]);
493
494 if(writeTime_)
495 fprintf(fp_," %4.4g ", update->dt*update->ntimestep);
496
497 fprintf(fp_," %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g %4.4g",
498 2.*radius[iPart]/force->cg(atom->type[iPart]),
499 x[iPart][0],x[iPart][1],x[iPart][2],
500 v[iPart][0],v[iPart][1],v[iPart][2]);
501
502 if(fix_orientation_)
503 {
504 double **orientation = NULL;
505 orientation = fix_orientation_->array_atom;
506 fprintf(fp_," %4.4g %4.4g %4.4g ",
507 orientation[iPart][0], orientation[iPart][1], orientation[iPart][2]);
508 }
509 if (fixColFound)
510 fprintf(fp_," %4.0g ", fix_color->vector_atom[iPart]);
511 fprintf(fp_,"\n");
512 fflush(fp_);
513 }
514 }
515
516 if(!delete_atoms_) //only set if not marked for deletion
517 counter[iPart] = once_ ? 2. : 1.;
518
519 }
520 else if(dot <= 0.) // dot <= 0
521 {
522 counter[iPart] = 0.;
523
524 }
525 }
526 }
527
528 MPI_Sum_Scalar(mass_this,world);
529 MPI_Sum_Scalar(nparticles_this,world);
530 if(fix_property_) MPI_Sum_Scalar(property_this,world);
531
532 mass_ += mass_this;
533 nparticles_ += nparticles_this;
534 property_sum_ += property_this;
535
536 }
537
538 /* ----------------------------------------------------------------------
539 perform particle deletion of marked particles
540 done before exchange, borders, reneighbor
541 so that ghost atoms and neighbor lists will be correct
542 ------------------------------------------------------------------------- */
543
pre_exchange()544 void FixMassflowMesh::pre_exchange()
545 {
546 if (delete_atoms_)
547 {
548 double mass_deleted_this_ = 0.;
549 double nparticles_deleted_this = 0.;
550 int *atom_map_array = atom->get_map_array();
551
552 // delete particles
553
554 while (atom_tags_delete_.size() > 0)
555 {
556 int iPart = atom->map(atom_tags_delete_[0]);
557
558 mass_deleted_this_ += atom->rmass[iPart];
559 nparticles_deleted_this += (fix_volumeweight_ms_ ? fix_volumeweight_ms_->vector_atom[iPart] : 1.);
560
561 atom->avec->copy(atom->nlocal-1,iPart,1);
562
563 atom_map_array[atom->tag[atom->nlocal-1]] = iPart;
564
565 atom->nlocal--;
566
567 atom_tags_delete_.erase(atom_tags_delete_.begin());
568 }
569
570 atom_tags_delete_.clear();
571
572 MPI_Sum_Scalar(mass_deleted_this_,world);
573 MPI_Sum_Scalar(nparticles_deleted_this,world);
574
575 mass_deleted_ += mass_deleted_this_;
576 nparticles_deleted_ += nparticles_deleted_this;
577
578 if(nparticles_deleted_this)
579 {
580 if (atom->tag_enable) {
581 if (atom->map_style) {
582 atom->nghost = 0;
583 atom->map_init();
584 atom->map_set();
585 }
586 }
587 }
588 }
589 }
590
591 /* ----------------------------------------------------------------------
592 pack entire state of Fix into one write
593 ------------------------------------------------------------------------- */
594
write_restart(FILE * fp)595 void FixMassflowMesh::write_restart(FILE *fp)
596 {
597 int n = 0;
598 double list[6];
599 list[n++] = mass_;
600 list[n++] = t_count_;
601 list[n++] = mass_last_;
602 list[n++] = nparticles_last_;
603 list[n++] = mass_deleted_;
604 list[n++] = nparticles_deleted_;
605
606 if (comm->me == 0) {
607 int size = n * sizeof(double);
608 fwrite(&size,sizeof(int),1,fp);
609 fwrite(list,sizeof(double),n,fp);
610 }
611 }
612
613 /* ----------------------------------------------------------------------
614 use state info from restart file to restart the Fix
615 ------------------------------------------------------------------------- */
616
restart(char * buf)617 void FixMassflowMesh::restart(char *buf)
618 {
619 int n = 0;
620 double *list = (double *) buf;
621
622 mass_ = list[n++];
623 t_count_ = list[n++];
624 mass_last_ = list[n++];
625 nparticles_last_ = list[n++];
626 mass_deleted_ = list[n++];
627 nparticles_deleted_ = list[n++];
628 }
629
630 /* ----------------------------------------------------------------------
631 output mass
632 ------------------------------------------------------------------------- */
633
compute_vector(int index)634 double FixMassflowMesh::compute_vector(int index)
635 {
636 if(reset_t_count_)
637 {
638 delta_t_ = t_count_;
639 t_count_ = 0.;
640 reset_t_count_ = false;
641 }
642
643 if(index == 0)
644 return mass_;
645 if(index == 1)
646 return nparticles_;
647 if(index == 2)
648 return delta_t_ == 0. ? 0. : (mass_-mass_last_)/delta_t_;
649 if(index == 3)
650 return delta_t_ == 0. ? 0. : (nparticles_-nparticles_last_)/delta_t_;
651 if(index == 4)
652 return mass_deleted_;
653 if(index == 5)
654 return nparticles_deleted_;
655 if(index == 6 && fix_property_)
656 return property_sum_;
657
658 return 0.;
659 }
660
661 /* ----------------------------------------------------------------------
662 get reference point on face
663 calculate normalvec
664 ------------------------------------------------------------------------- */
665
setRefPoint()666 void FixMassflowMesh::setRefPoint()
667 {
668 fix_mesh_->triMesh()->node(0,0,pref_);
669 fix_mesh_->triMesh()->surfaceNorm(0,nvec_);
670 double dot = vectorDot3D(nvec_,sidevec_);
671
672 if(fabs(dot) < 1e-6 && !havePointAtOutlet_ )
673 error->fix_error(FLERR,this,"need to change 'vec_side', it is currently in or to close to the mesh plane \n"
674 "This error may be caused by a moving mesh command, since 'vec_side' is not moved with the mesh.");
675 else if(dot < 0.)
676 vectorScalarMult3D(nvec_,-1.);
677 }
678