Monday, November 12, 2012

Triangular region in LAMMPS

I tried to use Steve Plimpton's  code for triangular region style, however it doesn't seem to work with the recent version of LAMMPS. I did minor modifications and publish them here.

src/region_triangle.h
#ifdef REGION_CLASS
RegionStyle(triangle,RegTriangle)
#else

#ifndef LMP_REGION_TRIANGLE_H
#define LMP_REGION_TRIANGLE_H

#include "region.h"

namespace LAMMPS_NS {

class RegTriangle : public Region {
 public:
  RegTriangle(class LAMMPS *, int, char **);
  ~RegTriangle();
  int inside(double, double, double);
  int surface_interior(double *, double);
  int surface_exterior(double *, double);

 private:
  double x1,y1,x2,y2,x3,y3;
};

}

#endif
#endif


src/region_triangle.cpp
#include "stdlib.h"
#include "string.h"
#include "region_triangle.h"
#include "error.h"

using namespace LAMMPS_NS;

#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))

/* ---------------------------------------------------------------------- */

RegTriangle::RegTriangle(LAMMPS *lmp, int narg, char **arg) :
  Region(lmp, narg, arg)
{
  options(narg-8,&arg[8]);

  x1 = xscale*atof(arg[2]);
  y1 = yscale*atof(arg[3]);
  x2 = xscale*atof(arg[4]);
  y2 = yscale*atof(arg[5]);
  x3 = xscale*atof(arg[6]);
  y3 = yscale*atof(arg[7]);

  // extent of 2d triangle

  extent_xlo = MIN(x1,x2);
  extent_xlo = MIN(extent_xlo,x3);
  extent_xhi = MAX(x1,x2);
  extent_xhi = MAX(extent_xhi,x3);

  extent_ylo = MIN(y1,y2);
  extent_ylo = MIN(extent_ylo,y3);
  extent_yhi = MAX(y1,y2);
  extent_yhi = MAX(extent_yhi,y3);

  extent_zlo = -0.5;
  extent_zhi = 0.5;

  cmax = 1;
  contact = new Contact[cmax];
}

 RegTriangle::~RegTriangle()
{
  delete [] contact;
}
 
/* ---------------------------------------------------------------------- */
 
int RegTriangle::inside(double x, double y, double z)
{
  double side1 = (x-x1)*(y2-y1) - (y-y1)*(x2-x1);
  double side2 = (x-x2)*(y3-y2) - (y-y2)*(x3-x2);
  double side3 = (x-x3)*(y1-y3) - (y-y3)*(x1-x3);
 
  if (side1 > 0.0 && side2 > 0.0 && side3 > 0.0) return 1;
  return 0;
}

/* ------------------------------------------------------------------------ */
int RegTriangle::surface_interior(double *x, double cutoff)
{
  if (RegTriangle::inside(x[0], x[1], x[2])==1)
   {
    return 1;
   }
  return 0;
}

/* ------------------------------------------------------------------------ */

int RegTriangle::surface_exterior(double *x, double cutoff)
{
  if (RegTriangle::inside(x[0], x[1], x[2])==1)
   {
    return 0;
   }
  return 1;
}


Notes and compilation
Functions surface_interior() and surface_exterior() are written for compatibility reasons and don't have much sense in present form. They are needed for fix/wall, so you should modify them if you want to use that fix. If not (as myself), just don't care.

Compilation is straightforward: put these two files in src directory and recompile LAMMPS. No further modification of the code is needed.

Example of use
Hex lattice with triangular void:

dimension       2
boundary        f m p

atom_style      atomic
lattice         hex 0.93
region          box block 0 300 0 500 -0.25 0.25
create_box      1 box
create_atoms    1 box
mass            1 1.0
region          lower block INF 18 INF INF INF INF units lattice
region          upper block 282 INF INF INF INF INF units lattice
group           lower region lower
group           upper region upper
group           boundary union lower upper
group           mobile subtract all boundary

# void
region          notch triangle 150 0 153 150 156 0
delete_atoms    region notch

No comments:

Post a Comment

Note: Only a member of this blog may post a comment.