git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11952 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/fix.h

@@ -101,6 +101,7 @@ class Fix : protected Pointers {
 
   virtual int setmask() = 0;
 
+  virtual void post_constructor() {}
   virtual void init() {}
   virtual void init_list(int, class NeighList *) {}
   virtual void setup(int) {}

src/fix_adapt.cpp

@@ -17,11 +17,13 @@
 #include "fix_adapt.h"
 #include "atom.h"
 #include "update.h"
+#include "group.h"
 #include "modify.h"
 #include "force.h"
 #include "pair.h"
 #include "pair_hybrid.h"
 #include "kspace.h"
+#include "fix_store.h"
 #include "input.h"
 #include "variable.h"
 #include "math_const.h"
@@ -152,6 +154,8 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   for (int m = 0; m < nadapt; m++)
     if (adapt[m].which == PAIR)
       memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
+
+  id_fix_diam = id_fix_chg = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -167,6 +171,13 @@ FixAdapt::~FixAdapt()
     }
   }
   delete [] adapt;
+
+  // check nfix in case all fixes have already been deleted
+
+  if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
+  if (id_fix_chg && modify->nfix) modify->delete_fix(id_fix_chg);
+  delete [] id_fix_diam;
+  delete [] id_fix_chg;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -179,6 +190,76 @@ int FixAdapt::setmask()
   return mask;
 }
 
+/* ----------------------------------------------------------------------
+   if need to restore per-atom quantities, create new fix STORE styles
+------------------------------------------------------------------------- */
+
+void FixAdapt::post_constructor()
+{
+  if (!resetflag) return;
+  if (!diamflag && !chgflag) return;
+
+  // new id = fix-ID + FIX_STORE_ATTRIBUTE
+  // new fix group = group for this fix
+
+  id_fix_diam = NULL;
+  id_fix_chg = NULL;
+
+  char **newarg = new char*[5];
+  newarg[1] = group->names[igroup];
+  newarg[2] = (char *) "STORE";
+  newarg[3] = (char *) "1";
+  newarg[4] = (char *) "1";
+
+  if (diamflag) {
+    int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
+    id_fix_diam = new char[n];
+    strcpy(id_fix_diam,id);
+    strcat(id_fix_diam,"_FIX_STORE_DIAM");
+    newarg[0] = id_fix_diam;
+    modify->add_fix(5,newarg);
+    fix_diam = (FixStore *) modify->fix[modify->nfix-1];
+
+    if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
+    else {
+      double *vec = fix_diam->vstore;
+      double *radius = atom->radius;
+      int *mask = atom->mask;
+      int nlocal = atom->nlocal;
+
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) vec[i] = radius[i];
+        else vec[i] = 0.0;
+      }
+    }
+  }
+
+  if (chgflag) {
+    int n = strlen(id) + strlen("_FIX_STORE_CHG") + 1;
+    id_fix_chg = new char[n];
+    strcpy(id_fix_chg,id);
+    strcat(id_fix_chg,"_FIX_STORE_CHG");
+    newarg[0] = id_fix_chg;
+    modify->add_fix(5,newarg);
+    fix_chg = (FixStore *) modify->fix[modify->nfix-1];
+
+    if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
+    else {
+      double *vec = fix_chg->vstore;
+      double *q = atom->q;
+      int *mask = atom->mask;
+      int nlocal = atom->nlocal;
+
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) vec[i] = q[i];
+        else vec[i] = 0.0;
+      }
+    }
+  }
+
+  delete [] newarg;
+}
+
 /* ---------------------------------------------------------------------- */
 
 void FixAdapt::init()
@@ -258,6 +339,19 @@ void FixAdapt::init()
           ad->array_orig[i][j] = ad->array[i][j];
     }
   }
+
+  // fixes that store initial per-atom values
+  
+  if (id_fix_diam) {
+    int ifix = modify->find_fix(id_fix_diam);
+    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
+    fix_diam = (FixStore *) modify->fix[ifix];
+  }
+  if (id_fix_chg) {
+    int ifix = modify->find_fix(id_fix_chg);
+    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
+    fix_chg = (FixStore *) modify->fix[ifix];
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -321,9 +415,11 @@ void FixAdapt::change_settings()
     } else if (ad->which == KSPACE) {
       *kspace_scale = value;
 
+    // set per atom values, also make changes for ghost atoms
+
     } else if (ad->which == ATOM) {
 
-      // set radius from diameter
+      // reset radius from diameter
       // also scale rmass to new value
 
       if (ad->aparam == DIAMETER) {
@@ -335,13 +431,14 @@ void FixAdapt::change_settings()
         double *rmass = atom->rmass;
         int *mask = atom->mask;
         int nlocal = atom->nlocal;
+        int nall = nlocal + atom->nghost;
 
         if (mflag == 0) {
-          for (i = 0; i < nlocal; i++)
+          for (i = 0; i < nall; i++)
             if (mask[i] & groupbit)
               radius[i] = 0.5*value;
         } else {
-          for (i = 0; i < nlocal; i++)
+          for (i = 0; i < nall; i++)
             if (mask[i] & groupbit) {
               density = rmass[i] / (4.0*MY_PI/3.0 *
                                     radius[i]*radius[i]*radius[i]);
@@ -354,7 +451,9 @@ void FixAdapt::change_settings()
         double *q = atom->q; 
 	int *mask = atom->mask;
 	int nlocal = atom->nlocal;
-	for (i = 0; i < nlocal; i++)
+        int nall = nlocal + atom->nghost;
+
+	for (i = 0; i < nall; i++)
 	  if (mask[i] & groupbit) q[i] = value; 
       }
     }
@@ -389,7 +488,34 @@ void FixAdapt::restore_settings()
       *kspace_scale = 1.0;
 
     } else if (ad->which == ATOM) {
+      if (diamflag) {
+        int mflag = 0;
+        if (atom->rmass_flag) mflag = 1;
+        double density;
 
+        double *vec = fix_diam->vstore;
+        double *radius = atom->radius;
+        double *rmass = atom->rmass;
+        int *mask = atom->mask;
+        int nlocal = atom->nlocal;
+
+        for (int i = 0; i < nlocal; i++)
+          if (mask[i] & groupbit) {
+            density = rmass[i] / (4.0*MY_PI/3.0 *
+                                  radius[i]*radius[i]*radius[i]);
+            radius[i] = vec[i];
+            rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density;
+          }
+      }
+      if (chgflag) {
+        double *vec = fix_chg->vstore;
+        double *q = atom->q;
+        int *mask = atom->mask;
+        int nlocal = atom->nlocal;
+
+        for (int i = 0; i < nlocal; i++)
+          if (mask[i] & groupbit) q[i] = vec[i];
+      }
     }
   }
 

src/fix_adapt.h

@@ -32,6 +32,7 @@ class FixAdapt : public Fix {
   FixAdapt(class LAMMPS *, int, char **);
   ~FixAdapt();
   int setmask();
+  void post_constructor();
   void init();
   void setup_pre_force(int);
   void pre_force(int);
@@ -40,6 +41,8 @@ class FixAdapt : public Fix {
  private:
   int nadapt,resetflag,scaleflag;
   int anypair;
+  char *id_fix_diam,*id_fix_chg;
+  class FixStore *fix_diam,*fix_chg;
 
   struct Adapt {
     int which,ivar;

src/modify.cpp

@@ -730,11 +730,6 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
 
   if (fix[ifix] == NULL) error->all(FLERR,"Invalid fix style");
 
-  // set fix mask values and increment nfix (if new)
-
-  fmask[ifix] = fix[ifix]->setmask();
-  if (newflag) nfix++;
-
   // check if Fix is in restart_global list
   // if yes, pass state info to the Fix so it can reset itself
 
@@ -766,6 +761,16 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
         if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
       }
     }
+
+  // increment nfix (if new)
+  // set fix mask values
+  // post_construct() allows new fix to create other fixes
+  // nfix increment comes first so that recursive call to add_fix within
+  //   post_constructor() will see updated nfix
+
+  if (newflag) nfix++;
+  fmask[ifix] = fix[ifix]->setmask();
+  fix[ifix]->post_constructor();
 }
 
 /* ----------------------------------------------------------------------