From 5a79c0b4ba72d69c5af5e4199d5e23a98f0a6fde Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 4 Dec 2008 16:14:49 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2281
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/delete_atoms.cpp | 33 +++++++++++++++++++++++++++++++++
 src/delete_atoms.h   |  1 +
 2 files changed, 34 insertions(+)

diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 03c6cf6f20..a567a6cfd2 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -24,6 +24,7 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "random_mars.h"
 #include "memory.h"
 #include "error.h"
 
@@ -56,6 +57,7 @@ void DeleteAtoms::command(int narg, char **arg)
   if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
   else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
   else if (strcmp(arg[0],"overlap") == 0) delete_overlap(narg,arg);
+  else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
   else error->all("Illegal delete_atoms command");
 
   // delete local atoms flagged in dlist
@@ -301,3 +303,34 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   create porosity by deleting atoms in a specified region 
+------------------------------------------------------------------------- */
+
+void DeleteAtoms::delete_porosity(int narg, char **arg)
+{
+  if (narg != 4) error->all("Illegal delete_atoms command");
+
+  int iregion = domain->find_region(arg[1]);
+  if (iregion == -1) error->all("Could not find delete_atoms region ID");
+
+  double porosity_fraction = atof(arg[2]);
+  int seed = atoi(arg[3]);
+
+  RanMars *random = new RanMars(lmp,seed + comm->me);
+
+  // allocate and initialize deletion list
+ 
+  int nlocal = atom->nlocal;
+  dlist = (int *) memory->smalloc(nlocal*sizeof(int),"delete_atoms:dlist");
+  for (int i = 0; i < nlocal; i++) dlist[i] = 0;
+
+  double **x = atom->x;
+
+  for (int i = 0; i < nlocal; i++)
+    if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
+      if (random->uniform() <= porosity_fraction) dlist[i] = 1;
+}
+
+
diff --git a/src/delete_atoms.h b/src/delete_atoms.h
index d929e4456d..2ccaca35f9 100644
--- a/src/delete_atoms.h
+++ b/src/delete_atoms.h
@@ -29,6 +29,7 @@ class DeleteAtoms : protected Pointers {
   void delete_group(int, char **);
   void delete_region(int, char **);
   void delete_overlap(int, char **);
+  void delete_porosity(int, char **);
 };
 
 }