/* * fix_printVolume.cpp * * Created on: Apr 18, 2012 * Author: kirill lykov */ #include "fix_printVolume.h" #include "stdlib.h" #include "string.h" #include "update.h" #include "input.h" #include "modify.h" #include "variable.h" #include "error.h" #include "atom.h" #include "neighbor.h" #include "domain.h" #include #include "math_const.h" #include "math_extra.h" using namespace LAMMPS_NS; namespace { void computeMassCenterTriangle(const double* xi1, const double* xi2, const double* xi3, double* massCenter) { for (int i = 0; i < 3; ++i) { massCenter[i] = (xi1[i] + xi2[i] + xi3[i]) / 3.; } } double computeTotalVolume(double** x, int* image, int** anglelist, int nanglelist, Domain*& domain, MPI_Comm comm) { double volume = 0.0; for (int n = 0; n < nanglelist; n++) { int i1 = anglelist[n][0]; int i2 = anglelist[n][1]; int i3 = anglelist[n][2]; // 1st bond - 21 double x21[3]; MathExtra::sub3(x[i2], x[i1], x21); // 3nd bond - 31 double x31[3]; MathExtra::sub3(x[i3], x[i1], x31); // normal to the triangle double normal[3]; MathExtra::cross3(x21, x31, normal); double massCenter[3]; computeMassCenterTriangle(x[i1], x[i2], x[i3], massCenter); domain->minimum_image(massCenter); double mcOnNorm = MathExtra::dot3(massCenter, normal); assert(mcOnNorm > 0.0); //if it is < 0 then mesh is not oriented properly volume += mcOnNorm / 6.; } double totalVolume; MPI_Allreduce(&volume, &totalVolume, 1, MPI_DOUBLE, MPI_SUM, comm); assert(totalVolume > 0.); //check that mesh is positively oriented return totalVolume; } } /* ---------------------------------------------------------------------- */ FixPrintVolume::FixPrintVolume(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 4) error->all(FLERR,"Illegal fix print command"); MPI_Comm_rank(world,&me); nevery = atoi(arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal fix print command"); } /* ---------------------------------------------------------------------- */ FixPrintVolume::~FixPrintVolume() { } /* ---------------------------------------------------------------------- */ int FixPrintVolume::setmask() { int mask = 0; mask |= END_OF_STEP; return mask; } /* ---------------------------------------------------------------------- */ void FixPrintVolume::end_of_step() { double** x = atom->x; int** anglelist = neighbor->anglelist; int nanglelist = neighbor->nanglelist; double res = computeTotalVolume(x, atom->image, anglelist, nanglelist, domain, world); if (me == 0) printf("Volume: %g", res); }