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

This commit is contained in:
sjplimp
2015-08-04 20:56:27 +00:00
parent 2f7344a02a
commit 318aefa19c
73 changed files with 15111 additions and 1 deletions

64
src/USER-SMD/Install.sh Normal file
View File

@ -0,0 +1,64 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*smd[^ \t]* //g' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I..\/..\/lib\/smd |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L..\/..\/lib\/smd$(LIBOBJDIR) |' ../Makefile.package
#sed -i -e 's|^PKG_LIB =[ \t]*|&-lsmd |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(smd_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(smd_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(smd_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*smd.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
# sed -i -e '4 i \
#include ..\/..\/lib\/smd\/Makefile.lammps
#' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*smd[^ \t]* //g' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*smd.*$/d' ../Makefile.package.settings
fi
fi

7
src/USER-SMD/README Normal file
View File

@ -0,0 +1,7 @@
Smooth Mach Dynamics -- Smooth Particle Hydrodynamics for solids and fluids
The person who created this package is Georg Ganzenmuller at the
Fraunhofer-Institute for High-Speed Dynamics, Ernst Mach Institute in
Germany (georg.ganzenmueller at emi.fhg.de). Contact him directly if
you have questions.

File diff suppressed because it is too large Load Diff

125
src/USER-SMD/atom_vec_smd.h Normal file
View File

@ -0,0 +1,125 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef ATOM_CLASS
AtomStyle(smd,AtomVecSMD)
#else
#ifndef LMP_ATOM_VEC_SMD_H
#define LMP_ATOM_VEC_SMD_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecSMD : public AtomVec {
public:
AtomVecSMD(class LAMMPS *);
~AtomVecSMD() {}
void init();
void grow(int);
void grow_reset();
void copy(int, int, int);
void force_clear(int, size_t);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
int unpack_comm_hybrid(int, int, double *);
int pack_reverse(int, int, double *);
int pack_reverse_hybrid(int, int, double *);
void unpack_reverse(int, int *, double *);
int unpack_reverse_hybrid(int, int *, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_vel(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
void unpack_border(int, int, double *);
void unpack_border_vel(int, int, double *);
int unpack_border_hybrid(int, int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **);
void data_vel(int, char **);
int data_vel_hybrid(int, char **);
void pack_data(double **);
int pack_data_hybrid(int, double *);
void write_data(FILE *, int, double **);
int write_data_hybrid(FILE *, double *);
void pack_vel(double **);
int pack_vel_hybrid(int, double *);
void write_vel(FILE *, int, double **);
int write_vel_hybrid(FILE *, double *);
bigint memory_usage();
private:
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f;
double *radius,*rmass;
int *molecule;
double *vfrac,**x0,*contact_radius, **smd_data_9, *e, *de, **vest;
double **tlsph_stress;
double *eff_plastic_strain;
double *damage;
double *eff_plastic_strain_rate;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Per-processor system is too big
The number of owned atoms plus ghost atoms on a single
processor must fit in 32-bit integer.
E: Invalid atom type in Atoms section of data file
Atom types must range from 1 to specified # of types.
E: Invalid radius in Atoms section of data file
Radius must be >= 0.0.
E: Invalid density in Atoms section of data file
Density value cannot be <= 0.0.
*/

View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_contact_radius.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDContactRadius::ComputeSMDContactRadius(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute smd/contact_radius command");
if (atom->contact_radius_flag != 1) error->all(FLERR,"compute smd/contact_radius command requires atom_style with contact_radius (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
contact_radius_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDContactRadius::~ComputeSMDContactRadius()
{
memory->sfree(contact_radius_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDContactRadius::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"smd/contact_radius") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute smd/contact_radius");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDContactRadius::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(contact_radius_vector);
nmax = atom->nmax;
contact_radius_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:contact_radius_vector");
vector_atom = contact_radius_vector;
}
double *contact_radius = atom->contact_radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
contact_radius_vector[i] = contact_radius[i];
}
else {
contact_radius_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDContactRadius::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/contact_radius,ComputeSMDContactRadius)
#else
#ifndef LMP_COMPUTE_SMD_CONTACT_RADIUS_H
#define LMP_COMPUTE_SMD_CONTACT_RADIUS_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDContactRadius : public Compute {
public:
ComputeSMDContactRadius(class LAMMPS *, int, char **);
~ComputeSMDContactRadius();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *contact_radius_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_damage.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDDamage::ComputeSMDDamage(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute smd/damage command");
if (atom->damage_flag != 1) error->all(FLERR,"compute smd/damage command requires atom_style with damage (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
damage_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDDamage::~ComputeSMDDamage()
{
memory->sfree(damage_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDDamage::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"smd/damage") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute smd/damage");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDDamage::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(damage_vector);
nmax = atom->nmax;
damage_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:damage_vector");
vector_atom = damage_vector;
}
double *damage = atom->damage;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
damage_vector[i] = damage[i];
}
else {
damage_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDDamage::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/damage,ComputeSMDDamage)
#else
#ifndef LMP_COMPUTE_SMD_DAMAGE_H
#define LMP_COMPUTE_SMD_DAMAGE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDDamage : public Compute {
public:
ComputeSMDDamage(class LAMMPS *, int, char **);
~ComputeSMDDamage();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *damage_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,112 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_hourglass_error.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDHourglassError::ComputeSMDHourglassError(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/hourglass_error command");
if (atom->smd_flag != 1)
error->all(FLERR, "compute smd/hourglass_error command requires atom_style with hourglass_error (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
hourglass_error_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDHourglassError::~ComputeSMDHourglassError() {
memory->sfree(hourglass_error_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDHourglassError::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/hourglass_error") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/hourglass_error");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDHourglassError::compute_peratom() {
invoked_peratom = update->ntimestep;
// grow output Vector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(hourglass_error_vector);
nmax = atom->nmax;
hourglass_error_vector = (double *) memory->smalloc(nmax * sizeof(double), "atom:hourglass_error_vector");
vector_atom = hourglass_error_vector;
}
int itmp = 0;
double *hourglass_error = (double *) force->pair->extract("smd/tlsph/hourglass_error_ptr", itmp);
if (hourglass_error == NULL) {
error->all(FLERR, "compute smd/hourglass_error failed to access hourglass_error array");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
hourglass_error_vector[i] = hourglass_error[i];
} else {
hourglass_error_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDHourglassError::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/hourglass_error,ComputeSMDHourglassError)
#else
#ifndef LMP_COMPUTE_SMD_HOURGLASS_ERROR_H
#define LMP_COMPUTE_SMD_HOURGLASS_ERROR_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDHourglassError : public Compute {
public:
ComputeSMDHourglassError(class LAMMPS *, int, char **);
~ComputeSMDHourglassError();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *hourglass_error_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_internal_energy.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDInternalEnergy::ComputeSMDInternalEnergy(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute smd/internal_energy command");
if (atom->e_flag != 1) error->all(FLERR,"compute smd/internal_energy command requires atom_style with internal_energy (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
internal_energy_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDInternalEnergy::~ComputeSMDInternalEnergy()
{
memory->sfree(internal_energy_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDInternalEnergy::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"smd/internal_energy") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute smd/internal_energy");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDInternalEnergy::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(internal_energy_vector);
nmax = atom->nmax;
internal_energy_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:internal_energy_vector");
vector_atom = internal_energy_vector;
}
double *e = atom->e;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
internal_energy_vector[i] = e[i];
}
else {
internal_energy_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDInternalEnergy::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/internal_energy,ComputeSMDInternalEnergy)
#else
#ifndef LMP_COMPUTE_SMD_INTERNAL_ENERGY_H
#define LMP_COMPUTE_SMD_INTERNAL_ENERGY_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDInternalEnergy : public Compute {
public:
ComputeSMDInternalEnergy(class LAMMPS *, int, char **);
~ComputeSMDInternalEnergy();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *internal_energy_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_plastic_strain.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDPlasticStrain::ComputeSMDPlasticStrain(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute smd/plastic_strain command");
if (atom->eff_plastic_strain_flag != 1) error->all(FLERR,"compute smd/plastic_strain command requires atom_style with plastic_strain (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
plastic_strain_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDPlasticStrain::~ComputeSMDPlasticStrain()
{
memory->sfree(plastic_strain_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDPlasticStrain::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"smd/plastic_strain") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute smd/plastic_strain");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDPlasticStrain::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(plastic_strain_vector);
nmax = atom->nmax;
plastic_strain_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:plastic_strain_vector");
vector_atom = plastic_strain_vector;
}
double *plastic_strain = atom->eff_plastic_strain;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
plastic_strain_vector[i] = plastic_strain[i];
}
else {
plastic_strain_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDPlasticStrain::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/plastic_strain,ComputeSMDPlasticStrain)
#else
#ifndef LMP_COMPUTE_SMD_PLASTIC_STRAIN_H
#define LMP_COMPUTE_SMD_PLASTIC_STRAIN_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDPlasticStrain : public Compute {
public:
ComputeSMDPlasticStrain(class LAMMPS *, int, char **);
~ComputeSMDPlasticStrain();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *plastic_strain_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_plastic_strain_rate.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDPlasticStrainRate::ComputeSMDPlasticStrainRate(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute smd/plastic_strain command");
if (atom->eff_plastic_strain_rate_flag != 1) error->all(FLERR,"compute smd/plastic_strain_rate command requires atom_style with plastic_strain_rate (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
plastic_strain_rate_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDPlasticStrainRate::~ComputeSMDPlasticStrainRate()
{
memory->sfree(plastic_strain_rate_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDPlasticStrainRate::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"smd/plastic_strain_rate") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute smd/plastic_strain_rate");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDPlasticStrainRate::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(plastic_strain_rate_vector);
nmax = atom->nmax;
plastic_strain_rate_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:plastic_strain_rate_vector");
vector_atom = plastic_strain_rate_vector;
}
double *plastic_strain_rate = atom->eff_plastic_strain_rate;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
plastic_strain_rate_vector[i] = plastic_strain_rate[i];
}
else {
plastic_strain_rate_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDPlasticStrainRate::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/plastic_strain_rate,ComputeSMDPlasticStrainRate)
#else
#ifndef LMP_COMPUTE_SMD_PLASTIC_STRAIN_RATE_H
#define LMP_COMPUTE_SMD_PLASTIC_STRAIN_RATE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDPlasticStrainRate : public Compute {
public:
ComputeSMDPlasticStrainRate(class LAMMPS *, int, char **);
~ComputeSMDPlasticStrainRate();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *plastic_strain_rate_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,107 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_rho.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDRho::ComputeSMDRho(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/rho command");
if (atom->vfrac_flag != 1)
error->all(FLERR, "compute smd/rho command requires atom_style with volume (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
rhoVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDRho::~ComputeSMDRho() {
memory->sfree(rhoVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDRho::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/rho") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/rho");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDRho::compute_peratom() {
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(rhoVector);
nmax = atom->nmax;
rhoVector = (double *) memory->smalloc(nmax * sizeof(double), "atom:rhoVector");
vector_atom = rhoVector;
}
double *vfrac = atom->vfrac;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
rhoVector[i] = rmass[i] / vfrac[i];
} else {
rhoVector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDRho::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/rho,ComputeSMDRho)
#else
#ifndef LMP_COMPUTE_SMD_RHO_H
#define LMP_COMPUTE_SMD_RHO_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDRho : public Compute {
public:
ComputeSMDRho(class LAMMPS *, int, char **);
~ComputeSMDRho();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *rhoVector;
};
}
#endif
#endif

View File

@ -0,0 +1,138 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_defgrad.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <iostream>
#include <stdio.h>
#include "stdlib.h"
#include "string.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHDefgrad::ComputeSMDTLSPHDefgrad(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_defgrad command");
peratom_flag = 1;
size_peratom_cols = 10;
nmax = 0;
defgradVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHDefgrad::~ComputeSMDTLSPHDefgrad() {
memory->sfree(defgradVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHDefgrad::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_defgrad") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_defgrad");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHDefgrad::compute_peratom() {
double **defgrad = atom->smd_data_9;
Matrix3d F;
invoked_peratom = update->ntimestep;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(defgradVector);
nmax = atom->nmax;
memory->create(defgradVector, nmax, size_peratom_cols, "defgradVector");
array_atom = defgradVector;
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
F(0, 0) = defgrad[i][0];
F(0, 1) = defgrad[i][1];
F(0, 2) = defgrad[i][2];
F(1, 0) = defgrad[i][3];
F(1, 1) = defgrad[i][4];
F(1, 2) = defgrad[i][5];
F(2, 0) = defgrad[i][6];
F(2, 1) = defgrad[i][7];
F(2, 2) = defgrad[i][8];
defgradVector[i][0] = F(0, 0);
defgradVector[i][1] = F(0, 1);
defgradVector[i][2] = F(0, 2);
defgradVector[i][3] = F(1, 0);
defgradVector[i][4] = F(1, 1);
defgradVector[i][5] = F(1, 2);
defgradVector[i][6] = F(2, 0);
defgradVector[i][7] = F(2, 1);
defgradVector[i][8] = F(2, 2);
defgradVector[i][9] = F.determinant();
} else {
defgradVector[i][0] = 1.0;
defgradVector[i][1] = 0.0;
defgradVector[i][2] = 0.0;
defgradVector[i][3] = 0.0;
defgradVector[i][4] = 1.0;
defgradVector[i][5] = 0.0;
defgradVector[i][6] = 0.0;
defgradVector[i][7] = 0.0;
defgradVector[i][8] = 1.0;
defgradVector[i][9] = 1.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTLSPHDefgrad::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_defgrad,ComputeSMDTLSPHDefgrad)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_DEFGRAD_H
#define LMP_COMPUTE_SMD_TLSPH_DEFGRAD_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDTLSPHDefgrad : public Compute {
public:
ComputeSMDTLSPHDefgrad(class LAMMPS *, int, char **);
~ComputeSMDTLSPHDefgrad();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double **defgradVector;
};
}
#endif
#endif

View File

@ -0,0 +1,115 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_dt.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTlsphDt::ComputeSMDTlsphDt(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_dt command");
if (atom->contact_radius_flag != 1)
error->all(FLERR,
"compute smd/tlsph_dt command requires atom_style with contact_radius (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
dt_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTlsphDt::~ComputeSMDTlsphDt() {
memory->sfree(dt_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTlsphDt::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_dt") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_dt");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTlsphDt::compute_peratom() {
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(dt_vector);
nmax = atom->nmax;
dt_vector = (double *) memory->smalloc(nmax * sizeof(double),
"atom:tlsph_dt_vector");
vector_atom = dt_vector;
}
int itmp = 0;
double *particle_dt = (double *) force->pair->extract("smd/tlsph/particle_dt_ptr",
itmp);
if (particle_dt == NULL) {
error->all(FLERR,
"compute smd/tlsph_dt failed to access particle_dt array");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dt_vector[i] = particle_dt[i];
} else {
dt_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTlsphDt::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_dt,ComputeSMDTlsphDt)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_DT_H
#define LMP_COMPUTE_SMD_TLSPH_DT_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDTlsphDt : public Compute {
public:
ComputeSMDTlsphDt(class LAMMPS *, int, char **);
~ComputeSMDTlsphDt();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *dt_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,107 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_num_neighs.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHNumNeighs::ComputeSMDTLSPHNumNeighs(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_num_neighs command");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
numNeighsRefConfigOutput = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHNumNeighs::~ComputeSMDTLSPHNumNeighs() {
memory->destroy(numNeighsRefConfigOutput);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHNumNeighs::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_num_neighs") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_num_neighs");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHNumNeighs::compute_peratom() {
invoked_peratom = update->ntimestep;
if (atom->nlocal > nmax) {
memory->destroy(numNeighsRefConfigOutput);
nmax = atom->nmax;
memory->create(numNeighsRefConfigOutput, nmax, "tlsph/num_neighs:numNeighsRefConfigOutput");
vector_atom = numNeighsRefConfigOutput;
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
int itmp = 0;
int *numNeighsRefConfig = (int *) force->pair->extract("smd/tlsph/numNeighsRefConfig_ptr", itmp);
if (numNeighsRefConfig == NULL) {
error->all(FLERR, "compute smd/tlsph_num_neighs failed to access numNeighsRefConfig array");
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
numNeighsRefConfigOutput[i] = numNeighsRefConfig[i];
} else {
numNeighsRefConfigOutput[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTLSPHNumNeighs::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,69 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_num_neighs,ComputeSMDTLSPHNumNeighs)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_NUM_NEIGHS_H
#define LMP_COMPUTE_SMD_TLSPH_NUM_NEIGHS_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDTLSPHNumNeighs : public Compute {
public:
ComputeSMDTLSPHNumNeighs(class LAMMPS *, int, char **);
~ComputeSMDTLSPHNumNeighs();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *numNeighsRefConfigOutput;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
W: More than one compute ke/atom
It is not efficient to use compute ke/atom more than once.
*/

View File

@ -0,0 +1,137 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_shape.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <iostream>
#include <stdio.h>
#include "stdlib.h"
#include "string.h"
#include <Eigen/Eigen>
#include <Eigen/Geometry>
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSmdTlsphShape::ComputeSmdTlsphShape(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_strain command");
peratom_flag = 1;
size_peratom_cols = 7;
nmax = 0;
strainVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSmdTlsphShape::~ComputeSmdTlsphShape() {
memory->sfree(strainVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSmdTlsphShape::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_strain") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_strain");
}
/* ---------------------------------------------------------------------- */
void ComputeSmdTlsphShape::compute_peratom() {
double *contact_radius = atom->contact_radius;
invoked_peratom = update->ntimestep;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(strainVector);
nmax = atom->nmax;
memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
array_atom = strainVector;
}
int itmp = 0;
Matrix3d *R = (Matrix3d *) force->pair->extract("smd/tlsph/rotation_ptr", itmp);
if (R == NULL) {
error->all(FLERR, "compute smd/tlsph_shape failed to access rotation array");
}
Matrix3d *F = (Matrix3d *) force->pair->extract("smd/tlsph/Fincr_ptr", itmp);
if (F == NULL) {
error->all(FLERR, "compute smd/tlsph_shape failed to access deformation gradient array");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
Matrix3d E, eye;
eye.setIdentity();
Quaterniond q;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
E = 0.5 * (F[i].transpose() * F[i] - eye); // Green-Lagrange strain
strainVector[i][0] = contact_radius[i] * (1.0 + E(0, 0));
strainVector[i][1] = contact_radius[i] * (1.0 + E(1, 1));
strainVector[i][2] = contact_radius[i] * (1.0 + E(2, 2));
q = R[i]; // convert pure rotation matrix to quaternion
strainVector[i][3] = q.w();
strainVector[i][4] = q.x();
strainVector[i][5] = q.y();
strainVector[i][6] = q.z();
} else {
for (int j = 0; j < size_peratom_cols; j++) {
strainVector[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSmdTlsphShape::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_shape,ComputeSmdTlsphShape)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_SHAPE_H
#define LMP_COMPUTE_SMD_TLSPH_SHAPE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSmdTlsphShape : public Compute {
public:
ComputeSmdTlsphShape(class LAMMPS *, int, char **);
~ComputeSmdTlsphShape();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double **strainVector;
};
}
#endif
#endif

View File

@ -0,0 +1,145 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_strain.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <iostream>
#include <stdio.h>
#include "stdlib.h"
#include "string.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHstrain::ComputeSMDTLSPHstrain(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_strain command");
peratom_flag = 1;
size_peratom_cols = 6;
nmax = 0;
strainVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHstrain::~ComputeSMDTLSPHstrain() {
memory->sfree(strainVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHstrain::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_strain") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_strain");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHstrain::compute_peratom() {
double **defgrad0 = atom->smd_data_9;
invoked_peratom = update->ntimestep;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(strainVector);
nmax = atom->nmax;
memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
array_atom = strainVector;
}
// copy data to output array
int itmp = 0;
Matrix3d *Fincr = (Matrix3d *) force->pair->extract("smd/tlsph/Fincr_ptr", itmp);
if (Fincr == NULL) {
error->all(FLERR, "compute smd/tlsph_strain failed to access Fincr array");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
Matrix3d E, eye, Ftotal, F0;
eye.setIdentity();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
// old deformation gradient
F0(0, 0) = defgrad0[i][0];
F0(0, 1) = defgrad0[i][1];
F0(0, 2) = defgrad0[i][2];
F0(1, 0) = defgrad0[i][3];
F0(1, 1) = defgrad0[i][4];
F0(1, 2) = defgrad0[i][5];
F0(2, 0) = defgrad0[i][6];
F0(2, 1) = defgrad0[i][7];
F0(2, 2) = defgrad0[i][8];
// compute current total deformation gradient
Ftotal = F0 * Fincr[i]; // this is the total deformation gradient: reference deformation times incremental deformation
E = 0.5 * (Ftotal.transpose() * Ftotal - eye); // Green-Lagrange strain
strainVector[i][0] = E(0, 0);
strainVector[i][1] = E(1, 1);
strainVector[i][2] = E(2, 2);
strainVector[i][3] = E(0, 1);
strainVector[i][4] = E(0, 2);
strainVector[i][5] = E(1, 2);
} else {
for (int j = 0; j < size_peratom_cols; j++) {
strainVector[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTLSPHstrain::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_strain,ComputeSMDTLSPHstrain)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_STRAIN_H
#define LMP_COMPUTE_SMD_TLSPH_STRAIN_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDTLSPHstrain : public Compute {
public:
ComputeSMDTLSPHstrain(class LAMMPS *, int, char **);
~ComputeSMDTLSPHstrain();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double **strainVector;
};
}
#endif
#endif

View File

@ -0,0 +1,126 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_strain_rate.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStrainRate::ComputeSMDTLSPHStrainRate(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/ulsph_strain_rate command");
peratom_flag = 1;
size_peratom_cols = 6;
nmax = 0;
strain_rate_array = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStrainRate::~ComputeSMDTLSPHStrainRate() {
memory->sfree(strain_rate_array);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHStrainRate::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/ulsph_strain_rate") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/ulsph_strain_rate");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHStrainRate::compute_peratom() {
invoked_peratom = update->ntimestep;
int *mol = atom->molecule;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(strain_rate_array);
nmax = atom->nmax;
memory->create(strain_rate_array, nmax, size_peratom_cols, "stresstensorVector");
array_atom = strain_rate_array;
}
int ulsph_flag = 0;
int tlsph_flag = 0;
int itmp = 0;
Matrix3d *D = (Matrix3d *) force->pair->extract("smd/tlsph/strain_rate_ptr", itmp);
if (D == NULL) {
error->all(FLERR,
"compute smd/tlsph_strain_rate could not access strain rate. Are the matching pair styles present?");
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
strain_rate_array[i][0] = D[i](0, 0); // xx
strain_rate_array[i][1] = D[i](1, 1); // yy
strain_rate_array[i][2] = D[i](2, 2); // zz
strain_rate_array[i][3] = D[i](0, 1); // xy
strain_rate_array[i][4] = D[i](0, 2); // xz
strain_rate_array[i][5] = D[i](1, 2); // yz
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTLSPHStrainRate::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDTLSPHStrainRate::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_strain_rate,ComputeSMDTLSPHStrainRate)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_STRAIN_RATE_H
#define LMP_COMPUTE_SMD_TLSPH_STRAIN_RATE_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class ComputeSMDTLSPHStrainRate : public Compute {
public:
ComputeSMDTLSPHStrainRate(class LAMMPS *, int, char **);
~ComputeSMDTLSPHStrainRate();
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;
double **strain_rate_array;
};
}
#endif
#endif

View File

@ -0,0 +1,131 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_tlsph_stress.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStress::ComputeSMDTLSPHStress(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_stress command");
peratom_flag = 1;
size_peratom_cols = 7;
nmax = 0;
stress_array = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTLSPHStress::~ComputeSMDTLSPHStress() {
memory->sfree(stress_array);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHStress::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_stress") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_stress");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTLSPHStress::compute_peratom() {
invoked_peratom = update->ntimestep;
Matrix3d stress_deviator;
double von_mises_stress;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(stress_array);
nmax = atom->nmax;
memory->create(stress_array, nmax, size_peratom_cols, "stresstensorVector");
array_atom = stress_array;
}
int itmp = 0;
Matrix3d *T = (Matrix3d *) force->pair->extract("smd/tlsph/stressTensor_ptr", itmp);
if (T == NULL) {
error->all(FLERR, "compute smd/tlsph_stress could not access stress tensors. Are the matching pair styles present?");
}
int nlocal = atom->nlocal;
int *mask = atom->mask;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
stress_deviator = Deviator(T[i]);
von_mises_stress = sqrt(3. / 2.) * stress_deviator.norm();
stress_array[i][0] = T[i](0, 0); // xx
stress_array[i][1] = T[i](1, 1); // yy
stress_array[i][2] = T[i](2, 2); // zz
stress_array[i][3] = T[i](0, 1); // xy
stress_array[i][4] = T[i](0, 2); // xz
stress_array[i][5] = T[i](1, 2); // yz
stress_array[i][6] = von_mises_stress;
} else {
for (int j = 0; j < size_peratom_cols; j++) {
stress_array[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTLSPHStress::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDTLSPHStress::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/tlsph_stress, ComputeSMDTLSPHStress)
#else
#ifndef LMP_COMPUTE_SMD_TLSPH_STRESS_H
#define LMP_COMPUTE_SMD_TLSPH_STRESS_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class ComputeSMDTLSPHStress : public Compute {
public:
ComputeSMDTLSPHStress(class LAMMPS *, int, char **);
~ComputeSMDTLSPHStress();
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;
double **stress_array;
};
}
#endif
#endif

View File

@ -0,0 +1,130 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_triangle_mesh_vertices.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <iostream>
#include <stdio.h>
#include "stdlib.h"
#include "string.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDTriangleVertices::ComputeSMDTriangleVertices(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/triangle_vertices command");
peratom_flag = 1;
size_peratom_cols = 9;
nmax = 0;
outputVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDTriangleVertices::~ComputeSMDTriangleVertices() {
memory->sfree(outputVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTriangleVertices::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/triangle_vertices") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/triangle_vertices");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDTriangleVertices::compute_peratom() {
double **smd_data_9 = atom->smd_data_9;
tagint *mol = atom->molecule;
invoked_peratom = update->ntimestep;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(outputVector);
nmax = atom->nmax;
memory->create(outputVector, nmax, size_peratom_cols, "defgradVector");
array_atom = outputVector;
}
/*
* triangle vertices are stored using the smd_data_9 array ...
* this is a hack but ok for now as I do not have to create additional storage space
* all triangle particles have molecule id >= 65535
*/
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && (mol[i] >= 65535) ){
outputVector[i][0] = smd_data_9[i][0];
outputVector[i][1] = smd_data_9[i][1];
outputVector[i][2] = smd_data_9[i][2];
outputVector[i][3] = smd_data_9[i][3];
outputVector[i][4] = smd_data_9[i][4];
outputVector[i][5] = smd_data_9[i][5];
outputVector[i][6] = smd_data_9[i][6];
outputVector[i][7] = smd_data_9[i][7];
outputVector[i][8] = smd_data_9[i][8];
} else {
for (int j = 0; j < size_peratom_cols; j++) {
outputVector[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDTriangleVertices::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/triangle_vertices,ComputeSMDTriangleVertices)
#else
#ifndef LMP_COMPUTE_SMD_TRIANGLE_VERTICES_H
#define LMP_COMPUTE_SMD_TRIANGLE_VERTICES_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDTriangleVertices : public Compute {
public:
ComputeSMDTriangleVertices(class LAMMPS *, int, char **);
~ComputeSMDTriangleVertices();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double **outputVector;
};
}
#endif
#endif

View File

@ -0,0 +1,115 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_ulsph_effm.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMD_Ulsph_Effm::ComputeSMD_Ulsph_Effm(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/ulsph_effm command");
if (atom->contact_radius_flag != 1)
error->all(FLERR,
"compute smd/ulsph_effm command requires atom_style with contact_radius (e.g. smd)");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
dt_vector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMD_Ulsph_Effm::~ComputeSMD_Ulsph_Effm() {
memory->sfree(dt_vector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMD_Ulsph_Effm::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/ulsph_effm") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/ulsph_effm");
}
/* ---------------------------------------------------------------------- */
void ComputeSMD_Ulsph_Effm::compute_peratom() {
invoked_peratom = update->ntimestep;
// grow rhoVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(dt_vector);
nmax = atom->nmax;
dt_vector = (double *) memory->smalloc(nmax * sizeof(double),
"atom:tlsph_dt_vector");
vector_atom = dt_vector;
}
int itmp = 0;
double *particle_dt = (double *) force->pair->extract("smd/ulsph/effective_modulus_ptr",
itmp);
if (particle_dt == NULL) {
error->all(FLERR,
"compute smd/ulsph_effm failed to access particle_dt array");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dt_vector[i] = particle_dt[i];
} else {
dt_vector[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMD_Ulsph_Effm::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/ulsph_effm,ComputeSMD_Ulsph_Effm)
#else
#ifndef LMP_COMPUTE_SMD_ULSPH_EFFM_H
#define LMP_COMPUTE_SMD_ULSPH_EFFM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMD_Ulsph_Effm : public Compute {
public:
ComputeSMD_Ulsph_Effm(class LAMMPS *, int, char **);
~ComputeSMD_Ulsph_Effm();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *dt_vector;
};
}
#endif
#endif

View File

@ -0,0 +1,107 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_ulsph_num_neighs.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHNumNeighs::ComputeSMDULSPHNumNeighs(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/ulsph_num_neighs command");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
numNeighsOutput = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHNumNeighs::~ComputeSMDULSPHNumNeighs() {
memory->destroy(numNeighsOutput);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHNumNeighs::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/ulsph_num_neighs") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/ulsph_num_neighs");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHNumNeighs::compute_peratom() {
invoked_peratom = update->ntimestep;
if (atom->nlocal > nmax) {
memory->destroy(numNeighsOutput);
nmax = atom->nmax;
memory->create(numNeighsOutput, nmax, "ulsph/num_neighs:numNeighsRefConfigOutput");
vector_atom = numNeighsOutput;
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
int itmp = 0;
int *numNeighs = (int *) force->pair->extract("smd/ulsph/numNeighs_ptr", itmp);
if (numNeighs == NULL) {
error->all(FLERR, "compute smd/ulsph_num_neighs failed to access numNeighs array");
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
numNeighsOutput[i] = numNeighs[i];
} else {
numNeighsOutput[i] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDULSPHNumNeighs::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,69 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/ulsph_num_neighs,ComputeSMDULSPHNumNeighs)
#else
#ifndef LMP_COMPUTE_SMD_ULSPH_NUM_NEIGHS_H
#define LMP_COMPUTE_SMD_ULSPH_NUM_NEIGHS_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDULSPHNumNeighs : public Compute {
public:
ComputeSMDULSPHNumNeighs(class LAMMPS *, int, char **);
~ComputeSMDULSPHNumNeighs();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *numNeighsOutput;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
W: More than one compute ke/atom
It is not efficient to use compute ke/atom more than once.
*/

View File

@ -0,0 +1,120 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_ulsph_strain.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <iostream>
#include <stdio.h>
#include "stdlib.h"
#include "string.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace std;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHstrain::ComputeSMDULSPHstrain(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/tlsph_strain command");
peratom_flag = 1;
size_peratom_cols = 6;
nmax = 0;
strainVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHstrain::~ComputeSMDULSPHstrain() {
memory->sfree(strainVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHstrain::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/tlsph_strain") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/tlsph_strain");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHstrain::compute_peratom() {
double **atom_data9 = atom->smd_data_9; // ULSPH strain is stored in the first 6 entries of this data field
invoked_peratom = update->ntimestep;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(strainVector);
nmax = atom->nmax;
memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
array_atom = strainVector;
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
strainVector[i][0] = atom_data9[i][0];
strainVector[i][1] = atom_data9[i][1];
strainVector[i][2] = atom_data9[i][2];
strainVector[i][3] = atom_data9[i][3];
strainVector[i][4] = atom_data9[i][4];
strainVector[i][5] = atom_data9[i][5];
} else {
for (int j = 0; j < size_peratom_cols; j++) {
strainVector[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDULSPHstrain::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/ulsph_strain,ComputeSMDULSPHstrain)
#else
#ifndef LMP_COMPUTE_SMD_ULSPH_STRAIN_H
#define LMP_COMPUTE_SMD_ULSPH_STRAIN_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDULSPHstrain : public Compute {
public:
ComputeSMDULSPHstrain(class LAMMPS *, int, char **);
~ComputeSMDULSPHstrain();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double **strainVector;
};
}
#endif
#endif

View File

@ -0,0 +1,132 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_ulsph_strain_rate.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStrainRate::ComputeSMDULSPHStrainRate(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/ulsph_strain_rate command");
peratom_flag = 1;
size_peratom_cols = 6;
nmax = 0;
strain_rate_array = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStrainRate::~ComputeSMDULSPHStrainRate() {
memory->sfree(strain_rate_array);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHStrainRate::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/ulsph_strain_rate") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/ulsph_strain_rate");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHStrainRate::compute_peratom() {
invoked_peratom = update->ntimestep;
int *mask = atom->mask;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(strain_rate_array);
nmax = atom->nmax;
memory->create(strain_rate_array, nmax, size_peratom_cols, "stresstensorVector");
array_atom = strain_rate_array;
}
int itmp = 0;
Matrix3d *L = (Matrix3d *) force->pair->extract("smd/ulsph/velocityGradient_ptr", itmp);
if (L == NULL) {
error->all(FLERR,
"compute smd/ulsph_strain_rate could not access any velocity gradients. Are the matching pair styles present?");
}
int nlocal = atom->nlocal;
Matrix3d D;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
D = 0.5 * (L[i] + L[i].transpose());
strain_rate_array[i][0] = D(0, 0); // xx
strain_rate_array[i][1] = D(1, 1); // yy
strain_rate_array[i][2] = D(2, 2); // zz
strain_rate_array[i][3] = D(0, 1); // xy
strain_rate_array[i][4] = D(0, 2); // xz
strain_rate_array[i][5] = D(1, 2); // yz
} else {
strain_rate_array[i][0] = 0.0;
strain_rate_array[i][1] = 0.0;
strain_rate_array[i][2] = 0.0;
strain_rate_array[i][3] = 0.0;
strain_rate_array[i][4] = 0.0;
strain_rate_array[i][5] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDULSPHStrainRate::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDULSPHStrainRate::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/ulsph_strain_rate,ComputeSMDULSPHStrainRate)
#else
#ifndef LMP_COMPUTE_SMD_ULSPH_STRAIN_RATE_H
#define LMP_COMPUTE_SMD_ULSPH_STRAIN_RATE_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class ComputeSMDULSPHStrainRate : public Compute {
public:
ComputeSMDULSPHStrainRate(class LAMMPS *, int, char **);
~ComputeSMDULSPHStrainRate();
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;
double **strain_rate_array;
};
}
#endif
#endif

View File

@ -0,0 +1,133 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_ulsph_stress.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStress::ComputeSMDULSPHStress(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/ulsph_stress command");
peratom_flag = 1;
size_peratom_cols = 7;
nmax = 0;
stress_array = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDULSPHStress::~ComputeSMDULSPHStress() {
memory->sfree(stress_array);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHStress::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/ulsph_stress") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/ulsph_stress");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDULSPHStress::compute_peratom() {
invoked_peratom = update->ntimestep;
int *mask = atom->mask;
Matrix3d stress_deviator;
double von_mises_stress;
// grow vector array if necessary
if (atom->nlocal > nmax) {
memory->destroy(stress_array);
nmax = atom->nmax;
memory->create(stress_array, nmax, size_peratom_cols, "stresstensorVector");
array_atom = stress_array;
}
int itmp = 0;
Matrix3d *T = (Matrix3d *) force->pair->extract("smd/ulsph/stressTensor_ptr", itmp);
if (T == NULL) {
error->all(FLERR, "compute smd/ulsph_stress could not access stress tensors. Are the matching pair styles present?");
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
stress_deviator = Deviator(T[i]);
von_mises_stress = sqrt(3. / 2.) * stress_deviator.norm();
stress_array[i][0] = T[i](0, 0); // xx
stress_array[i][1] = T[i](1, 1); // yy
stress_array[i][2] = T[i](2, 2); // zz
stress_array[i][3] = T[i](0, 1); // xy
stress_array[i][4] = T[i](0, 2); // xz
stress_array[i][5] = T[i](1, 2); // yz
stress_array[i][6] = von_mises_stress;
} else {
for (int j = 0; j < size_peratom_cols; j++) {
stress_array[i][j] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDULSPHStress::memory_usage() {
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}
/*
* deviator of a tensor
*/
Matrix3d ComputeSMDULSPHStress::Deviator(Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}

View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/ulsph_stress, ComputeSMDULSPHStress)
#else
#ifndef LMP_COMPUTE_SMD_ULSPH_STRESS_H
#define LMP_COMPUTE_SMD_ULSPH_STRESS_H
#include "compute.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class ComputeSMDULSPHStress : public Compute {
public:
ComputeSMDULSPHStress(class LAMMPS *, int, char **);
~ComputeSMDULSPHStress();
void init();
void compute_peratom();
double memory_usage();
Matrix3d Deviator(Matrix3d);
private:
int nmax;
double **stress_array;
};
}
#endif
#endif

View File

@ -0,0 +1,130 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_smd_vol.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSMDVol::ComputeSMDVol(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) {
if (narg != 3)
error->all(FLERR, "Illegal compute smd/volume command");
if (atom->vfrac_flag != 1)
error->all(FLERR, "compute smd/volume command requires atom_style with density (e.g. smd)");
scalar_flag = 1;
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
volVector = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeSMDVol::~ComputeSMDVol() {
memory->sfree(volVector);
}
/* ---------------------------------------------------------------------- */
void ComputeSMDVol::init() {
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style, "smd/volume") == 0)
count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR, "More than one compute smd/volume");
}
/* ---------------------------------------------------------------------- */
void ComputeSMDVol::compute_peratom() {
invoked_peratom = update->ntimestep;
// grow volVector array if necessary
if (atom->nlocal > nmax) {
memory->sfree(volVector);
nmax = atom->nmax;
volVector = (double *) memory->smalloc(nmax * sizeof(double), "atom:volVector");
vector_atom = volVector;
}
double *vfrac = atom->vfrac;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
volVector[i] = vfrac[i];
} else {
volVector[i] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
double ComputeSMDVol::compute_scalar() {
invoked_scalar = update->ntimestep;
double *vfrac = atom->vfrac;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double this_proc_sum_volumes = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
this_proc_sum_volumes += vfrac[i];
}
}
//printf("this_proc_sum_volumes = %g\n", this_proc_sum_volumes);
MPI_Allreduce(&this_proc_sum_volumes, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
//if (comm->me == 0) printf("global sum_volumes = %g\n", scalar);
return scalar;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeSMDVol::memory_usage() {
double bytes = nmax * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,56 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(smd/volume,ComputeSMDVol)
#else
#ifndef LMP_COMPUTE_SMD_VOL_H
#define LMP_COMPUTE_SMD_VOL_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSMDVol : public Compute {
public:
ComputeSMDVol(class LAMMPS *, int, char **);
~ComputeSMDVol();
void init();
void compute_peratom();
double compute_scalar();
double memory_usage();
private:
int nmax;
double *volVector;
};
}
#endif
#endif

View File

@ -0,0 +1,233 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_smd_adjust_dt.h"
#include "atom.h"
#include "update.h"
#include "integrate.h"
#include "domain.h"
#include "lattice.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "fix.h"
#include "output.h"
#include "dump.h"
#include "comm.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
FixSMDTlsphDtReset::FixSMDTlsphDtReset(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (narg != 4)
error->all(FLERR, "Illegal fix smd/adjust_dt command");
// set time_depend, else elapsed time accumulation can be messed up
time_depend = 1;
scalar_flag = 1;
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extscalar = 0;
extvector = 0;
restart_global = 1; // this fix stores global (i.e., not per-atom) info: elaspsed time
safety_factor = atof(arg[3]);
// initializations
t_elapsed = 0.0;
}
/* ---------------------------------------------------------------------- */
int FixSMDTlsphDtReset::setmask() {
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDTlsphDtReset::init() {
dt = update->dt;
}
/* ---------------------------------------------------------------------- */
void FixSMDTlsphDtReset::setup(int vflag) {
end_of_step();
}
/* ---------------------------------------------------------------------- */
void FixSMDTlsphDtReset::initial_integrate(int vflag) {
//printf("in adjust_dt: dt = %20.10f\n", update->dt);
t_elapsed += update->dt;
}
/* ---------------------------------------------------------------------- */
void FixSMDTlsphDtReset::end_of_step() {
double dtmin = BIG;
int itmp = 0;
/*
* extract minimum CFL timestep from TLSPH and ULSPH pair styles
*/
double *dtCFL_TLSPH = (double *) force->pair->extract("smd/tlsph/dtCFL_ptr", itmp);
double *dtCFL_ULSPH = (double *) force->pair->extract("smd/ulsph/dtCFL_ptr", itmp);
double *dt_TRI = (double *) force->pair->extract("smd/tri_surface/stable_time_increment_ptr", itmp);
double *dt_HERTZ = (double *) force->pair->extract("smd/hertz/stable_time_increment_ptr", itmp);
double *dt_PERI_IPMB = (double *) force->pair->extract("smd/peri_ipmb/stable_time_increment_ptr", itmp);
if ((dtCFL_TLSPH == NULL) && (dtCFL_ULSPH == NULL) && (dt_TRI == NULL) && (dt_HERTZ == NULL)
&& (dt_PERI_IPMB == NULL)) {
error->all(FLERR, "fix smd/adjust_dt failed to access a valid dtCFL");
}
if (dtCFL_TLSPH != NULL) {
dtmin = MIN(dtmin, *dtCFL_TLSPH);
}
if (dtCFL_ULSPH != NULL) {
dtmin = MIN(dtmin, *dtCFL_ULSPH);
}
if (dt_TRI != NULL) {
dtmin = MIN(dtmin, *dt_TRI);
}
if (dt_HERTZ != NULL) {
dtmin = MIN(dtmin, *dt_HERTZ);
}
if (dt_PERI_IPMB != NULL) {
dtmin = MIN(dtmin, *dt_PERI_IPMB);
}
// double **v = atom->v;
// double **f = atom->f;
// double *rmass = atom->rmass;
// double *radius = atom->radius;
// int *mask = atom->mask;
// int nlocal = atom->nlocal;
// double dtv, dtf, dtsq;
// double vsq, fsq, massinv, xmax;
// double delx, dely, delz, delr;
// for (int i = 0; i < nlocal; i++) {
// if (mask[i] & groupbit) {
// xmax = 0.005 * radius[i];
// massinv = 1.0 / rmass[i];
// vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
// fsq = f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2];
// dtv = dtf = BIG;
// if (vsq > 0.0)
// dtv = xmax / sqrt(vsq);
// if (fsq > 0.0)
// dtf = sqrt(2.0 * xmax / (sqrt(fsq) * massinv));
// dt = MIN(dtv, dtf);
// dtmin = MIN(dtmin, dt);
// dtsq = dt * dt;
// delx = dt * v[i][0] + 0.5 * dtsq * massinv * f[i][0];
// dely = dt * v[i][1] + 0.5 * dtsq * massinv * f[i][1];
// delz = dt * v[i][2] + 0.5 * dtsq * massinv * f[i][2];
// delr = sqrt(delx * delx + dely * dely + delz * delz);
// if (delr > xmax)
// dt *= xmax / delr;
// dtmin = MIN(dtmin, dt);
//
//// xmax = 0.05 * radius[i];
//// massinv = 1.0 / rmass[i];
//// fsq = f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2];
//// dtf = BIG;
//// if (fsq > 0.0)
//// dtf = sqrt(2.0 * xmax / (sqrt(fsq) * massinv));
//// dtmin = MIN(dtmin, dtf);
// }
// }
dtmin *= safety_factor; // apply safety factor
MPI_Allreduce(&dtmin, &dt, 1, MPI_DOUBLE, MPI_MIN, world);
if (update->ntimestep == 0) {
dt = 1.0e-16;
}
//printf("dtmin is now: %f, dt is now%f\n", dtmin, dt);
update->dt = dt;
if (force->pair)
force->pair->reset_dt();
for (int i = 0; i < modify->nfix; i++)
modify->fix[i]->reset_dt();
}
/* ---------------------------------------------------------------------- */
double FixSMDTlsphDtReset::compute_scalar() {
return t_elapsed;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixSMDTlsphDtReset::write_restart(FILE *fp) {
int n = 0;
double list[1];
list[n++] = t_elapsed;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size, sizeof(int), 1, fp);
fwrite(list, sizeof(double), n, fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixSMDTlsphDtReset::restart(char *buf) {
int n = 0;
double *list = (double *) buf;
t_elapsed = list[n++];
}

View File

@ -0,0 +1,80 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/adjust_dt,FixSMDTlsphDtReset)
#else
#ifndef LMP_FIX_TLSPH_DT_RESET_H
#define LMP_FIX_TLSPH_DT_RESET_H
#include "fix.h"
namespace LAMMPS_NS {
class FixSMDTlsphDtReset: public Fix {
public:
FixSMDTlsphDtReset(class LAMMPS *, int, char **);
~FixSMDTlsphDtReset() {
}
int setmask();
void init();
void setup(int);
void initial_integrate(int);
void end_of_step();
double compute_scalar();
void write_restart(FILE *);
void restart(char *);
private:
double safety_factor;
double dt, t_elapsed;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Use of fix dt/reset with undefined lattice
Must use lattice command with fix dt/reset command if units option is
set to lattice.
W: Dump dcd/xtc timestamp may be wrong with fix dt/reset
If the fix changes the timestep, the dump dcd file will not
reflect the change.
*/

View File

@ -0,0 +1,254 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_smd_integrate_tlsph.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "error.h"
#include "pair.h"
#include "neigh_list.h"
#include <Eigen/Eigen>
#include "domain.h"
#include "neighbor.h"
#include "comm.h"
#include "modify.h"
#include "stdio.h"
#include <iostream>
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace std;
/* ---------------------------------------------------------------------- */
FixSMDIntegrateTlsph::FixSMDIntegrateTlsph(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (narg < 3) {
printf("narg=%d\n", narg);
error->all(FLERR, "Illegal fix smd/integrate_tlsph command");
}
xsphFlag = false;
vlimit = -1.0;
int iarg = 3;
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix smd/integrate_tlsph is active for group: %s \n", arg[1]);
}
while (true) {
if (iarg >= narg) {
break;
}
if (strcmp(arg[iarg], "xsph") == 0) {
xsphFlag = true;
if (comm->me == 0) {
error->one(FLERR, "XSPH is currently not available");
printf("... will use XSPH time integration\n");
}
} else if (strcmp(arg[iarg], "limit_velocity") == 0) {
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected number following limit_velocity");
}
vlimit = force->numeric(FLERR, arg[iarg]);
if (comm->me == 0) {
printf("... will limit velocities to <= %g\n", vlimit);
}
} else {
char msg[128];
sprintf(msg, "Illegal keyword for smd/integrate_tlsph: %s\n", arg[iarg]);
error->all(FLERR, msg);
}
iarg++;
}
if (comm->me == 0) {
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
time_integrate = 1;
// set comm sizes needed by this fix
atom->add_callback(0);
}
/* ---------------------------------------------------------------------- */
int FixSMDIntegrateTlsph::setmask() {
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateTlsph::init() {
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
vlimitsq = vlimit * vlimit;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixSMDIntegrateTlsph::initial_integrate(int vflag) {
double dtfm, vsq, scale;
// update v and x of atoms in group
double **x = atom->x;
double **v = atom->v;
double **vest = atom->vest;
double **f = atom->f;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int itmp;
double vxsph_x, vxsph_y, vxsph_z;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
Vector3d *smoothVelDifference = (Vector3d *) force->pair->extract("smd/tlsph/smoothVel_ptr", itmp);
if (xsphFlag) {
if (smoothVelDifference == NULL) {
error->one(FLERR,
"fix smd/integrate_tlsph failed to access smoothVel array. Check if a pair style exist which calculates this quantity.");
}
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
// 1st part of Velocity_Verlet: push velocties 1/2 time increment ahead
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vlimit > 0.0) {
vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
if (vsq > vlimitsq) {
scale = sqrt(vlimitsq / vsq);
v[i][0] *= scale;
v[i][1] *= scale;
v[i][2] *= scale;
}
}
if (xsphFlag) {
// construct XSPH velocity
vxsph_x = v[i][0] + 0.5 * smoothVelDifference[i](0);
vxsph_y = v[i][1] + 0.5 * smoothVelDifference[i](1);
vxsph_z = v[i][2] + 0.5 * smoothVelDifference[i](2);
vest[i][0] = vxsph_x + dtfm * f[i][0];
vest[i][1] = vxsph_y + dtfm * f[i][1];
vest[i][2] = vxsph_z + dtfm * f[i][2];
x[i][0] += dtv * vxsph_x;
x[i][1] += dtv * vxsph_y;
x[i][2] += dtv * vxsph_z;
} else {
// extrapolate velocity from half- to full-step
vest[i][0] = v[i][0] + dtfm * f[i][0];
vest[i][1] = v[i][1] + dtfm * f[i][1];
vest[i][2] = v[i][2] + dtfm * f[i][2];
x[i][0] += dtv * v[i][0]; // 2nd part of Velocity-Verlet: push positions one full time increment ahead
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateTlsph::final_integrate() {
double dtfm, vsq, scale;
// update v of atoms in group
double **v = atom->v;
double **f = atom->f;
double *e = atom->e;
double *de = atom->de;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
int i;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0]; // 3rd part of Velocity-Verlet: push velocities another half time increment ahead
v[i][1] += dtfm * f[i][1]; // both positions and velocities are now defined at full time-steps.
v[i][2] += dtfm * f[i][2];
// limit velocity
if (vlimit > 0.0) {
vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
if (vsq > vlimitsq) {
scale = sqrt(vlimitsq / vsq);
v[i][0] *= scale;
v[i][1] *= scale;
v[i][2] *= scale;
}
}
e[i] += dtv * de[i];
}
}
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateTlsph::reset_dt() {
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
vlimitsq = vlimit * vlimit;
}
/* ---------------------------------------------------------------------- */

View File

@ -0,0 +1,73 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/integrate_tlsph,FixSMDIntegrateTlsph)
#else
#ifndef LMP_FIX_SMD_INTEGRATE_TLSPH_H
#define LMP_FIX_SMD_INTEGRATE_TLSPH_H
#include "fix.h"
namespace LAMMPS_NS {
class FixSMDIntegrateTlsph: public Fix {
friend class Neighbor;
friend class PairTlsph;
public:
FixSMDIntegrateTlsph(class LAMMPS *, int, char **);
virtual ~FixSMDIntegrateTlsph() {
}
int setmask();
virtual void init();
virtual void initial_integrate(int);
virtual void final_integrate();
virtual void reset_dt();
protected:
double dtv, dtf, vlimit, vlimitsq;
int mass_require;
bool xsphFlag;
class Pair *pair;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

View File

@ -0,0 +1,325 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdio.h"
#include "string.h"
#include "fix_smd_integrate_ulsph.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include "domain.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixSMDIntegrateUlsph::FixSMDIntegrateUlsph(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if ((atom->e_flag != 1) || (atom->vfrac_flag != 1))
error->all(FLERR, "fix smd/integrate_ulsph command requires atom_style with both energy and volume");
if (narg < 3)
error->all(FLERR, "Illegal number of arguments for fix smd/integrate_ulsph command");
adjust_radius_flag = false;
xsphFlag = false;
vlimit = -1.0;
int iarg = 3;
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix smd/integrate_ulsph is active for group: %s \n", arg[1]);
}
while (true) {
if (iarg >= narg) {
break;
}
if (strcmp(arg[iarg], "xsph") == 0) {
xsphFlag = true;
if (comm->me == 0) {
error->one(FLERR, "XSPH is currently not available");
printf("... will use XSPH time integration\n");
}
} else if (strcmp(arg[iarg], "adjust_radius") == 0) {
adjust_radius_flag = true;
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
}
adjust_radius_factor = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
}
min_nn = force->inumeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
}
max_nn = force->inumeric(FLERR, arg[iarg]);
if (comm->me == 0) {
printf("... will adjust smoothing length dynamically with factor %g to achieve %d to %d neighbors per particle.\n ",
adjust_radius_factor, min_nn, max_nn);
}
} else if (strcmp(arg[iarg], "limit_velocity") == 0) {
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected number following limit_velocity");
}
vlimit = force->numeric(FLERR, arg[iarg]);
if (comm->me == 0) {
printf("... will limit velocities to <= %g\n", vlimit);
}
} else {
char msg[128];
sprintf(msg, "Illegal keyword for smd/integrate_ulsph: %s\n", arg[iarg]);
error->all(FLERR, msg);
}
iarg++;
}
if (comm->me == 0) {
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
// set comm sizes needed by this fix
atom->add_callback(0);
time_integrate = 1;
}
/* ---------------------------------------------------------------------- */
int FixSMDIntegrateUlsph::setmask() {
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateUlsph::init() {
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
vlimitsq = vlimit * vlimit;
}
/* ----------------------------------------------------------------------
allow for both per-type and per-atom mass
------------------------------------------------------------------------- */
void FixSMDIntegrateUlsph::initial_integrate(int vflag) {
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **vest = atom->vest;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int i, itmp;
double dtfm, vsq, scale;
double vxsph_x, vxsph_y, vxsph_z;
//printf("initial_integrate at timestep %d\n", update->ntimestep);
/*
* get smoothed velocities from ULSPH pair style
*/
Vector3d *smoothVel = (Vector3d *) force->pair->extract("smd/ulsph/smoothVel_ptr", itmp);
if (xsphFlag) {
if (smoothVel == NULL) {
error->one(FLERR, "fix smd/integrate_ulsph failed to access smoothVel array");
}
}
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vlimit > 0.0) {
vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
if (vsq > vlimitsq) {
scale = sqrt(vlimitsq / vsq);
v[i][0] *= scale;
v[i][1] *= scale;
v[i][2] *= scale;
vest[i][0] = v[i][0];
vest[i][1] = v[i][1];
vest[i][2] = v[i][2];
}
}
if (xsphFlag) {
// construct XSPH velocity
vxsph_x = v[i][0] - 0.5 * smoothVel[i](0);
vxsph_y = v[i][1] - 0.5 * smoothVel[i](1);
vxsph_z = v[i][2] - 0.5 * smoothVel[i](2);
vest[i][0] = vxsph_x + dtfm * f[i][0];
vest[i][1] = vxsph_y + dtfm * f[i][1];
vest[i][2] = vxsph_z + dtfm * f[i][2];
x[i][0] += dtv * vxsph_x;
x[i][1] += dtv * vxsph_y;
x[i][2] += dtv * vxsph_z;
} else {
// extrapolate velocity from half- to full-step
vest[i][0] = v[i][0] + dtfm * f[i][0];
vest[i][1] = v[i][1] + dtfm * f[i][1];
vest[i][2] = v[i][2] + dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateUlsph::final_integrate() {
double **v = atom->v;
double **f = atom->f;
double *e = atom->e;
double *de = atom->de;
double *vfrac = atom->vfrac;
double *radius = atom->radius;
double *contact_radius = atom->contact_radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
double dtfm, vsq, scale;
double *rmass = atom->rmass;
double vol_increment;
Matrix3d D;
/*
* get current number of SPH neighbors from ULSPH pair style
*/
int itmp;
int *nn = (int *) force->pair->extract("smd/ulsph/numNeighs_ptr", itmp);
if (nn == NULL) {
error->one(FLERR, "fix smd/integrate_ulsph failed to accesss num_neighs array");
}
Matrix3d *L = (Matrix3d *) force->pair->extract("smd/ulsph/velocityGradient_ptr", itmp);
if (L == NULL) {
error->one(FLERR, "fix smd/integrate_ulsph failed to accesss velocityGradient array");
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vlimit > 0.0) {
vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
if (vsq > vlimitsq) {
scale = sqrt(vlimitsq / vsq);
v[i][0] *= scale;
v[i][1] *= scale;
v[i][2] *= scale;
}
}
e[i] += dtf * de[i];
if (adjust_radius_flag) {
if (nn[i] < min_nn) {
radius[i] *= adjust_radius_factor;
} else if (nn[i] > max_nn) {
radius[i] /= adjust_radius_factor;
}
radius[i] = MAX(radius[i], 1.25 * contact_radius[i]);
radius[i] = MIN(radius[i], 4.0 * contact_radius[i]);
}
D = 0.5 * (L[i] + L[i].transpose());
vol_increment = vfrac[i] * update->dt * D.trace(); // Jacobian of deformation
vfrac[i] += vol_increment;
}
}
}
/* ---------------------------------------------------------------------- */
void FixSMDIntegrateUlsph::reset_dt() {
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}

View File

@ -0,0 +1,64 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/integrate_ulsph,FixSMDIntegrateUlsph)
#else
#ifndef LMP_FIX_SMD_INTEGRATE_ULSPH_H
#define LMP_FIX_SMD_INTEGRATE_ULSPH_H
#include "fix.h"
namespace LAMMPS_NS {
class FixSMDIntegrateUlsph : public Fix {
public:
FixSMDIntegrateUlsph(class LAMMPS *, int, char **);
int setmask();
virtual void init();
virtual void initial_integrate(int);
virtual void final_integrate();
void reset_dt();
private:
class NeighList *list;
protected:
double dtv,dtf, vlimit, vlimitsq;;
int mass_require;
bool xsphFlag;
bool adjust_radius_flag;
double adjust_radius_factor;
int min_nn, max_nn; // number of SPH neighbors should lie within this interval
class Pair *pair;
};
}
#endif
#endif

View File

@ -0,0 +1,517 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdio.h"
#include "string.h"
#include "fix_smd_move_triangulated_surface.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include "domain.h"
#include <Eigen/Eigen>
#include "math_const.h"
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using namespace std;
/* ---------------------------------------------------------------------- */
FixSMDMoveTriSurf::FixSMDMoveTriSurf(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (atom->smd_flag != 1) {
error->all(FLERR, "fix fix smd/move_tri_surf command requires atom_style smd");
}
if (narg < 3)
error->all(FLERR, "Illegal number of arguments for fix fix smd/move_tri_surf command");
rotateFlag = linearFlag = wiggleFlag = false;
wiggle_direction = 1.0;
wiggle_max_travel = 0.0;
int iarg = 3;
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix fix smd/move_tri_surf is active for group: %s \n", arg[1]);
}
while (true) {
if (iarg >= narg) {
break;
}
if (strcmp(arg[iarg], "*LINEAR") == 0) {
linearFlag = true;
if (comm->me == 0) {
printf("... will move surface in a linear fashion\n");
}
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vx = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vy = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vz = force->numeric(FLERR, arg[iarg]);
} else if (strcmp(arg[iarg], "*WIGGLE") == 0) {
wiggleFlag = true;
if (comm->me == 0) {
printf("... will move surface in wiggle fashion\n");
}
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vx = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vy = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vz = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
wiggle_max_travel = force->numeric(FLERR, arg[iarg]);
} else if (strcmp(arg[iarg], "*ROTATE") == 0) {
rotateFlag = true;
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(0) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(1) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(2) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(0) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(1) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(2) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_period = force->numeric(FLERR, arg[iarg]);
/*
* construct rotation matrix
*/
u_cross(0, 0) = 0.0;
u_cross(0, 1) = -rotation_axis(2);
u_cross(0, 2) = rotation_axis(1);
u_cross(1, 0) = rotation_axis(2);
u_cross(1, 1) = 0.0;
u_cross(1, 2) = -rotation_axis(0);
u_cross(2, 0) = -rotation_axis(1);
u_cross(2, 1) = rotation_axis(0);
u_cross(2, 2) = 0.0;
uxu = rotation_axis * rotation_axis.transpose();
if (comm->me == 0) {
printf("will rotate with period %f\n", rotation_period);
}
} else {
char msg[128];
sprintf(msg, "Illegal keyword for fix smd/move_tri_surf: %s\n", arg[iarg]);
error->all(FLERR, msg);
}
iarg++;
}
if (comm->me == 0) {
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
// set comm sizes needed by this fix
comm_forward = 12;
//atom->add_callback(0);
//atom->add_callback(1);
time_integrate = 1;
}
/* ---------------------------------------------------------------------- */
FixSMDMoveTriSurf::~FixSMDMoveTriSurf()
{
// unregister callbacks to this fix from Atom class
//atom->delete_callback(id,0);
//atom->delete_callback(id,1);
}
/* ---------------------------------------------------------------------- */
int FixSMDMoveTriSurf::setmask() {
int mask = 0;
mask |= INITIAL_INTEGRATE;
//mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::init() {
dtv = update->dt;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixSMDMoveTriSurf::initial_integrate(int vflag) {
double **x = atom->x;
double **x0 = atom->x0;
double **v = atom->v;
double **vest = atom->vest;
double **smd_data_9 = atom->smd_data_9;
int *mol = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double phi;
int i;
Matrix3d eye, Rot;
eye.setIdentity();
Vector3d v1, v2, v3, n, point, rotated_point, R, vel;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
if (linearFlag) { // translate particles
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] = vx;
v[i][1] = vy;
v[i][2] = vz;
vest[i][0] = vx;
vest[i][1] = vy;
vest[i][2] = vz;
x[i][0] += dtv * vx;
x[i][1] += dtv * vy;
x[i][2] += dtv * vz;
/*
* if this is a triangle, move the vertices as well
*/
if (mol[i] >= 65535) {
smd_data_9[i][0] += dtv * vx;
smd_data_9[i][1] += dtv * vy;
smd_data_9[i][2] += dtv * vz;
smd_data_9[i][3] += dtv * vx;
smd_data_9[i][4] += dtv * vy;
smd_data_9[i][5] += dtv * vz;
smd_data_9[i][6] += dtv * vx;
smd_data_9[i][7] += dtv * vy;
smd_data_9[i][8] += dtv * vz;
}
}
}
}
if (wiggleFlag) { // wiggle particles forward and backward
wiggle_travel += sqrt(vx * vx + vy * vy + vz * vz ) * dtv;
double wiggle_vx = wiggle_direction * vx;
double wiggle_vy = wiggle_direction * vy;
double wiggle_vz = wiggle_direction * vz;
//printf("wiggle vz is %f, wiggle_max_travel is %f, dir=%f\n", wiggle_vz, wiggle_max_travel, wiggle_direction);
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] = wiggle_vx;
v[i][1] = wiggle_vy;
v[i][2] = wiggle_vz;
vest[i][0] = wiggle_vx;
vest[i][1] = wiggle_vy;
vest[i][2] = wiggle_vz;
x[i][0] += dtv * wiggle_vx;
x[i][1] += dtv * wiggle_vy;
x[i][2] += dtv * wiggle_vz;
/*
* if this is a triangle, move the vertices as well
*/
if (mol[i] >= 65535) {
smd_data_9[i][0] += dtv * wiggle_vx;
smd_data_9[i][1] += dtv * wiggle_vy;
smd_data_9[i][2] += dtv * wiggle_vz;
smd_data_9[i][3] += dtv * wiggle_vx;
smd_data_9[i][4] += dtv * wiggle_vy;
smd_data_9[i][5] += dtv * wiggle_vz;
smd_data_9[i][6] += dtv * wiggle_vx;
smd_data_9[i][7] += dtv * wiggle_vy;
smd_data_9[i][8] += dtv * wiggle_vz;
}
}
}
if (wiggle_travel >= wiggle_max_travel) {
wiggle_direction *= -1.0;
wiggle_travel = 0.0;
}
}
if (rotateFlag) { // rotate particles
Vector3d xnew, R_new, x_correct;
/*
* rotation angle and matrix form of Rodrigues' rotation formula
*/
phi = MY_2PI * dtv / rotation_period;
//printf("dt=%f, phi =%f, T=%f\n", dtv, phi, rotation_period);
Rot = cos(phi) * eye + sin(phi) * u_cross + (1.0 - cos(phi)) * uxu;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
/*
* generate vector R from origin to point which is to be rotated
*/
point << x[i][0], x[i][1], x[i][2];
R = point - origin;
/*
* rotate vector and shift away from origin
*/
rotated_point = Rot * R + origin;
/*
* determine velocity
*/
vel = (rotated_point - point) / update->dt;
/*
* assign new velocities and coordinates
*/
v[i][0] = vel(0);
v[i][1] = vel(1);
v[i][2] = vel(2);
vest[i][0] = vel(0);
vest[i][1] = vel(1);
vest[i][2] = vel(2);
x[i][0] = rotated_point(0);
x[i][1] = rotated_point(1);
x[i][2] = rotated_point(2);
/*
* if this is a triangle, rotate the vertices as well
*/
if (mol[i] >= 65535) {
v1 << smd_data_9[i][0], smd_data_9[i][1], smd_data_9[i][2];
R = v1 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v1) / update->dt;
smd_data_9[i][0] = rotated_point(0);
smd_data_9[i][1] = rotated_point(1);
smd_data_9[i][2] = rotated_point(2);
v1 = rotated_point;
v2 << smd_data_9[i][3], smd_data_9[i][4], smd_data_9[i][5];
R = v2 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v2) / update->dt;
smd_data_9[i][3] = rotated_point(0);
smd_data_9[i][4] = rotated_point(1);
smd_data_9[i][5] = rotated_point(2);
v2 = rotated_point;
v3 << smd_data_9[i][6], smd_data_9[i][7], smd_data_9[i][8];
R = v3 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v3) / update->dt;
smd_data_9[i][6] = rotated_point(0);
smd_data_9[i][7] = rotated_point(1);
smd_data_9[i][8] = rotated_point(2);
v3 = rotated_point;
// recalculate triangle normal
n = (v2 - v1).cross(v2 - v3);
n /= n.norm();
x0[i][0] = n(0);
x0[i][1] = n(1);
x0[i][2] = n(2);
}
}
}
}
// we changed smd_data_9, x0. perform communication to ghosts
comm->forward_comm_fix(this);
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::reset_dt() {
dtv = update->dt;
}
/* ---------------------------------------------------------------------- */
int FixSMDMoveTriSurf::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
int i, j, m;
double **x0 = atom->x0;
double **smd_data_9 = atom->smd_data_9;
//printf("in FixSMDIntegrateTlsph::pack_forward_comm\n");
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
buf[m++] = smd_data_9[j][0];
buf[m++] = smd_data_9[j][1];
buf[m++] = smd_data_9[j][2];
buf[m++] = smd_data_9[j][3];
buf[m++] = smd_data_9[j][4];
buf[m++] = smd_data_9[j][5];
buf[m++] = smd_data_9[j][6];
buf[m++] = smd_data_9[j][7];
buf[m++] = smd_data_9[j][8];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::unpack_forward_comm(int n, int first, double *buf) {
int i, m, last;
double **x0 = atom->x0;
double **smd_data_9 = atom->smd_data_9;
//printf("in FixSMDMoveTriSurf::unpack_forward_comm\n");
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
smd_data_9[i][0] = buf[m++];
smd_data_9[i][1] = buf[m++];
smd_data_9[i][2] = buf[m++];
smd_data_9[i][3] = buf[m++];
smd_data_9[i][4] = buf[m++];
smd_data_9[i][5] = buf[m++];
smd_data_9[i][6] = buf[m++];
smd_data_9[i][7] = buf[m++];
smd_data_9[i][8] = buf[m++];
}
}

View File

@ -0,0 +1,65 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/move_tri_surf,FixSMDMoveTriSurf)
#else
#ifndef LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H
#define LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H
#include "fix.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class FixSMDMoveTriSurf: public Fix {
public:
FixSMDMoveTriSurf(class LAMMPS *, int, char **);
~FixSMDMoveTriSurf();
int setmask();
virtual void init();
virtual void initial_integrate(int);
void reset_dt();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
protected:
double dtv;
bool linearFlag, wiggleFlag, rotateFlag;
double vx, vy, vz;
Vector3d rotation_axis, origin;
double rotation_period;
Matrix3d u_cross, uxu;
double wiggle_travel, wiggle_max_travel, wiggle_direction;
};
}
#endif
#endif

View File

@ -0,0 +1,375 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "fix_smd_setvel.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum {
NONE, CONSTANT, EQUAL, ATOM
};
/* ---------------------------------------------------------------------- */
FixSMDSetVel::FixSMDSetVel(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (narg < 6)
error->all(FLERR, "Illegal fix setvelocity command");
dynamic_group_allow = 1;
vector_flag = 1;
size_vector = 3;
global_freq = 1;
extvector = 1;
xstr = ystr = zstr = NULL;
if (strstr(arg[3], "v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];
strcpy(xstr, &arg[3][2]);
} else if (strcmp(arg[3], "NULL") == 0) {
xstyle = NONE;
} else {
xvalue = force->numeric(FLERR, arg[3]);
xstyle = CONSTANT;
}
if (strstr(arg[4], "v_") == arg[4]) {
int n = strlen(&arg[4][2]) + 1;
ystr = new char[n];
strcpy(ystr, &arg[4][2]);
} else if (strcmp(arg[4], "NULL") == 0) {
ystyle = NONE;
} else {
yvalue = force->numeric(FLERR, arg[4]);
ystyle = CONSTANT;
}
if (strstr(arg[5], "v_") == arg[5]) {
int n = strlen(&arg[5][2]) + 1;
zstr = new char[n];
strcpy(zstr, &arg[5][2]);
} else if (strcmp(arg[5], "NULL") == 0) {
zstyle = NONE;
} else {
zvalue = force->numeric(FLERR, arg[5]);
zstyle = CONSTANT;
}
// optional args
iregion = -1;
idregion = NULL;
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg], "region") == 0) {
if (iarg + 2 > narg)
error->all(FLERR, "Illegal fix setvelocity command");
iregion = domain->find_region(arg[iarg + 1]);
if (iregion == -1)
error->all(FLERR, "Region ID for fix setvelocity does not exist");
int n = strlen(arg[iarg + 1]) + 1;
idregion = new char[n];
strcpy(idregion, arg[iarg + 1]);
iarg += 2;
} else
error->all(FLERR, "Illegal fix setvelocity command");
}
force_flag = 0;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
maxatom = atom->nmax;
memory->create(sforce, maxatom, 3, "setvelocity:sforce");
}
/* ---------------------------------------------------------------------- */
FixSMDSetVel::~FixSMDSetVel() {
delete[] xstr;
delete[] ystr;
delete[] zstr;
delete[] idregion;
memory->destroy(sforce);
}
/* ---------------------------------------------------------------------- */
int FixSMDSetVel::setmask() {
int mask = 0;
//mask |= INITIAL_INTEGRATE;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDSetVel::init() {
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0)
error->all(FLERR, "Variable name for fix setvelocity does not exist");
if (input->variable->equalstyle(xvar))
xstyle = EQUAL;
else if (input->variable->atomstyle(xvar))
xstyle = ATOM;
else
error->all(FLERR, "Variable for fix setvelocity is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0)
error->all(FLERR, "Variable name for fix setvelocity does not exist");
if (input->variable->equalstyle(yvar))
ystyle = EQUAL;
else if (input->variable->atomstyle(yvar))
ystyle = ATOM;
else
error->all(FLERR, "Variable for fix setvelocity is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0)
error->all(FLERR, "Variable name for fix setvelocity does not exist");
if (input->variable->equalstyle(zvar))
zstyle = EQUAL;
else if (input->variable->atomstyle(zvar))
zstyle = ATOM;
else
error->all(FLERR, "Variable for fix setvelocity is invalid style");
}
// set index and check validity of region
if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR, "Region ID for fix setvelocity does not exist");
}
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else
varflag = CONSTANT;
if (strstr(update->integrate_style, "respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// cannot use non-zero forces for a minimization since no energy is integrated
// use fix addforce instead
int flag = 0;
if (update->whichflag == 2) {
if (xstyle == EQUAL || xstyle == ATOM)
flag = 1;
if (ystyle == EQUAL || ystyle == ATOM)
flag = 1;
if (zstyle == EQUAL || zstyle == ATOM)
flag = 1;
if (xstyle == CONSTANT && xvalue != 0.0)
flag = 1;
if (ystyle == CONSTANT && yvalue != 0.0)
flag = 1;
if (zstyle == CONSTANT && zvalue != 0.0)
flag = 1;
}
if (flag)
error->all(FLERR, "Cannot use non-zero forces in an energy minimization");
}
/* ---------------------------------------------------------------------- */
void FixSMDSetVel::setup(int vflag) {
if (strstr(update->integrate_style, "verlet"))
post_force(vflag);
else
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
((Respa *) update->integrate)->copy_flevel_f(ilevel);
post_force_respa(vflag, ilevel, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel);
}
}
/* ---------------------------------------------------------------------- */
void FixSMDSetVel::min_setup(int vflag) {
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
//void FixSMDSetVel::initial_integrate(int vflag) {
void FixSMDSetVel::post_force(int vflag) {
double **x = atom->x;
double **f = atom->f;
double **v = atom->v;
double **vest = atom->vest;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// update region if necessary
Region *region = NULL;
if (iregion >= 0) {
region = domain->regions[iregion];
region->prematch();
}
// reallocate sforce array if necessary
if (varflag == ATOM && nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(sforce);
memory->create(sforce, maxatom, 3, "setvelocity:sforce");
}
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
force_flag = 0;
if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2]))
continue;
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (xstyle) {
v[i][0] = xvalue;
vest[i][0] = xvalue;
f[i][0] = 0.0;
}
if (ystyle) {
v[i][1] = yvalue;
vest[i][1] = yvalue;
f[i][1] = 0.0;
}
if (zstyle) {
v[i][2] = zvalue;
vest[i][2] = zvalue;
f[i][2] = 0.0;
}
}
// variable force, wrap with clear/add
} else {
modify->clearstep_compute();
if (xstyle == EQUAL)
xvalue = input->variable->compute_equal(xvar);
else if (xstyle == ATOM)
input->variable->compute_atom(xvar, igroup, &sforce[0][0], 3, 0);
if (ystyle == EQUAL)
yvalue = input->variable->compute_equal(yvar);
else if (ystyle == ATOM)
input->variable->compute_atom(yvar, igroup, &sforce[0][1], 3, 0);
if (zstyle == EQUAL)
zvalue = input->variable->compute_equal(zvar);
else if (zstyle == ATOM)
input->variable->compute_atom(zvar, igroup, &sforce[0][2], 3, 0);
modify->addstep_compute(update->ntimestep + 1);
//printf("setting velocity at timestep %d\n", update->ntimestep);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2]))
continue;
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (xstyle == ATOM) {
vest[i][0] = v[i][0] = sforce[i][0];
f[i][0] = 0.0;
} else if (xstyle) {
vest[i][0] = v[i][0] = xvalue;
f[i][0] = 0.0;
}
if (ystyle == ATOM) {
vest[i][1] = v[i][1] = sforce[i][1];
f[i][1] = 0.0;
} else if (ystyle) {
vest[i][1] = v[i][1] = yvalue;
f[i][1] = 0.0;
}
if (zstyle == ATOM) {
vest[i][2] = v[i][2] = sforce[i][2];
f[i][2] = 0.0;
} else if (zstyle) {
vest[i][2] = v[i][2] = zvalue;
f[i][2] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
return components of total force on fix group before force was changed
------------------------------------------------------------------------- */
double FixSMDSetVel::compute_vector(int n) {
// only sum across procs one time
if (force_flag == 0) {
MPI_Allreduce(foriginal, foriginal_all, 3, MPI_DOUBLE, MPI_SUM, world);
force_flag = 1;
}
return foriginal_all[n];
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixSMDSetVel::memory_usage() {
double bytes = 0.0;
if (varflag == ATOM)
bytes = atom->nmax * 3 * sizeof(double);
return bytes;
}

View File

@ -0,0 +1,95 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/setvelocity,FixSMDSetVel)
#else
#ifndef LMP_FIX_SMD_SET_VELOCITY_H
#define LMP_FIX_SMD_SET_VELOCITY_H
#include "fix.h"
namespace LAMMPS_NS {
class FixSMDSetVel : public Fix {
public:
FixSMDSetVel(class LAMMPS *, int, char **);
~FixSMDSetVel();
int setmask();
void init();
void setup(int);
void min_setup(int);
//void initial_integrate(int);
void post_force(int);
double compute_vector(int);
double memory_usage();
private:
double xvalue,yvalue,zvalue;
int varflag,iregion;
char *xstr,*ystr,*zstr;
char *idregion;
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[3],foriginal_all[3];
int force_flag;
int nlevels_respa;
int maxatom;
double **sforce;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Region ID for fix setforce does not exist
Self-explanatory.
E: Variable name for fix setforce does not exist
Self-explanatory.
E: Variable for fix setforce is invalid style
Only equal-style variables can be used.
E: Cannot use non-zero forces in an energy minimization
Fix setforce cannot be used in this manner. Use fix addforce
instead.
*/

View File

@ -0,0 +1,613 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* This file is based on the FixShearHistory class.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "lattice.h"
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "fix_smd_tlsph_reference_configuration.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "pair.h"
#include "update.h"
#include "modify.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include <Eigen/Eigen>
#include "smd_kernels.h"
#include "smd_math.h"
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace SMD_Kernels;
using namespace std;
using namespace SMD_Math;
#define DELTA 16384
#define INSERT_PREDEFINED_CRACKS false
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::FixSMD_TLSPH_ReferenceConfiguration(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (atom->map_style == 0)
error->all(FLERR, "Pair tlsph with partner list requires an atom map, see atom_modify");
maxpartner = 1;
npartner = NULL;
partner = NULL;
wfd_list = NULL;
wf_list = NULL;
energy_per_bond = NULL;
degradation_ij = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// initialize npartner to 0 so neighbor list creation is OK the 1st time
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
npartner[i] = 0;
}
comm_forward = 14;
updateFlag = 1;
}
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::~FixSMD_TLSPH_ReferenceConfiguration() {
// unregister this fix so atom class doesn't invoke it any more
atom->delete_callback(id, 0);
// delete locally stored arrays
memory->destroy(npartner);
memory->destroy(partner);
memory->destroy(wfd_list);
memory->destroy(wf_list);
memory->destroy(degradation_ij);
memory->destroy(energy_per_bond);
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::setmask() {
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::init() {
if (atom->tag_enable == 0)
error->all(FLERR, "Pair style tlsph requires atoms have IDs");
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() {
//return;
//printf("in FixSMD_TLSPH_ReferenceConfiguration::pre_exchange()\n");
double **defgrad = atom->smd_data_9;
double *radius = atom->radius;
double *rho = atom->rho;
double *vfrac = atom->vfrac;
double **x = atom->x;
double **x0 = atom->x0;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int i, itmp;
int *mask = atom->mask;
if (igroup == atom->firstgroup) {
nlocal = atom->nfirst;
}
int *updateFlag_ptr = (int *) force->pair->extract("smd/tlsph/updateFlag_ptr", itmp);
if (updateFlag_ptr == NULL) {
error->one(FLERR,
"fix FixSMD_TLSPH_ReferenceConfiguration failed to access updateFlag pointer. Check if a pair style exist which calculates this quantity.");
}
int *nn = (int *) force->pair->extract("smd/tlsph/numNeighsRefConfig_ptr", itmp);
if (nn == NULL) {
error->all(FLERR, "FixSMDIntegrateTlsph::updateReferenceConfiguration() failed to access numNeighsRefConfig_ptr array");
}
// sum all update flag across processors
MPI_Allreduce(updateFlag_ptr, &updateFlag, 1, MPI_INT, MPI_MAX, world);
if (updateFlag > 0) {
if (comm->me == 0) {
printf("**** updating ref config at step: %ld\n", update->ntimestep);
}
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
// re-set x0 coordinates
x0[i][0] = x[i][0];
x0[i][1] = x[i][1];
x0[i][2] = x[i][2];
// re-set deformation gradient
defgrad[i][0] = 1.0;
defgrad[i][1] = 0.0;
defgrad[i][2] = 0.0;
defgrad[i][3] = 0.0;
defgrad[i][4] = 1.0;
defgrad[i][5] = 0.0;
defgrad[i][6] = 0.0;
defgrad[i][7] = 0.0;
defgrad[i][8] = 1.0;
/*
* Adjust particle volume as the reference configuration is changed.
* We safeguard against excessive deformations by limiting the adjustment range
* to the intervale J \in [0.9..1.1]
*/
vfrac[i] = rmass[i] / rho[i];
//
if (nn[i] < 15) {
radius[i] *= 1.2;
} // else //{
// radius[i] *= pow(J, 1.0 / domain->dimension);
//}
}
}
// update of reference config could have changed x0, vfrac, radius
// communicate these quantities now to ghosts: x0, vfrac, radius
comm->forward_comm_fix(this);
setup(0);
}
}
/* ----------------------------------------------------------------------
copy partner info from neighbor lists to atom arrays
so can be migrated or stored with atoms
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) {
int i, j, ii, jj, n, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
int itype, jtype;
double r, h, wf, wfd;
Vector3d dx;
if (updateFlag == 0)
return;
int nlocal = atom->nlocal;
nmax = atom->nmax;
grow_arrays(nmax);
// 1st loop over neighbor list
// calculate npartner for each owned atom
// nlocal_neigh = nlocal when neigh list was built, may be smaller than nlocal
double **x0 = atom->x;
double *radius = atom->radius;
int *type = atom->type;
int *mask = atom->mask;
tagint *tag = atom->tag;
NeighList *list = pair->list;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// zero npartner for all current atoms
for (i = 0; i < nlocal; i++)
npartner[i] = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (r <= h) {
npartner[i]++;
if (j < nlocal) {
npartner[j]++;
}
}
}
}
maxpartner = 0;
for (i = 0; i < nlocal; i++)
maxpartner = MAX(maxpartner, npartner[i]);
int maxall;
MPI_Allreduce(&maxpartner, &maxall, 1, MPI_INT, MPI_MAX, world);
maxpartner = maxall;
grow_arrays(nmax);
for (i = 0; i < nlocal; i++) {
npartner[i] = 0;
for (jj = 0; jj < maxpartner; jj++) {
wfd_list[i][jj] = 0.0;
wf_list[i][jj] = 0.0;
degradation_ij[i][jj] = 0.0;
energy_per_bond[i][jj] = 0.0;
}
}
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
jtype = type[j];
h = radius[i] + radius[j];
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
if (r < h) {
spiky_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
partner[i][npartner[i]] = tag[j];
wfd_list[i][npartner[i]] = wfd;
wf_list[i][npartner[i]] = wf;
npartner[i]++;
if (j < nlocal) {
partner[j][npartner[j]] = tag[i];
wfd_list[j][npartner[j]] = wfd;
wf_list[j][npartner[j]] = wf;
npartner[j]++;
}
}
}
}
// count number of particles for which this group is active
// bond statistics
if (update->ntimestep > -1) {
n = 0;
int count = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
n += npartner[i];
count += 1;
}
}
int nall, countall;
MPI_Allreduce(&n, &nall, 1, MPI_INT, MPI_SUM, world);
MPI_Allreduce(&count, &countall, 1, MPI_INT, MPI_SUM, world);
if (comm->me == 0) {
if (screen) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, "TLSPH neighbors:\n");
fprintf(screen, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(screen, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
if (logfile) {
fprintf(logfile, "\nTLSPH neighbors:\n");
fprintf(logfile, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(logfile, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
}
}
}
updateFlag = 0; // set update flag to zero after the update
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixSMD_TLSPH_ReferenceConfiguration::memory_usage() {
int nmax = atom->nmax;
int bytes = nmax * sizeof(int);
bytes += nmax * maxpartner * sizeof(tagint); // partner array
bytes += nmax * maxpartner * sizeof(float); // wf_list
bytes += nmax * maxpartner * sizeof(float); // wfd_list
bytes += nmax * maxpartner * sizeof(float); // damage_per_interaction array
bytes += nmax * sizeof(int); // npartner array
return bytes;
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::grow_arrays(int nmax) {
//printf("in FixSMD_TLSPH_ReferenceConfiguration::grow_arrays\n");
memory->grow(npartner, nmax, "tlsph_refconfig_neigh:npartner");
memory->grow(partner, nmax, maxpartner, "tlsph_refconfig_neigh:partner");
memory->grow(wfd_list, nmax, maxpartner, "tlsph_refconfig_neigh:wfd");
memory->grow(wf_list, nmax, maxpartner, "tlsph_refconfig_neigh:wf");
memory->grow(degradation_ij, nmax, maxpartner, "tlsph_refconfig_neigh:degradation_ij");
memory->grow(energy_per_bond, nmax, maxpartner, "tlsph_refconfig_neigh:damage_onset_strain");
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::copy_arrays(int i, int j, int delflag) {
npartner[j] = npartner[i];
for (int m = 0; m < npartner[j]; m++) {
partner[j][m] = partner[i][m];
wfd_list[j][m] = wfd_list[i][m];
wf_list[j][m] = wf_list[i][m];
degradation_ij[j][m] = degradation_ij[i][m];
energy_per_bond[j][m] = energy_per_bond[i][m];
}
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_exchange(int i, double *buf) {
// NOTE: how do I know comm buf is big enough if extreme # of touching neighs
// Comm::BUFEXTRA may need to be increased
//printf("pack_exchange ...\n");
int m = 0;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::unpack_exchange(int nlocal, double *buf) {
if (nlocal == nmax) {
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
nmax = nmax / DELTA * DELTA;
nmax += DELTA;
grow_arrays(nmax);
error->message(FLERR,
"in Fixtlsph_refconfigNeighGCG::unpack_exchange: local arrays too small for receiving partner information; growing arrays");
}
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
int m = 0;
npartner[nlocal] = static_cast<int>(buf[m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<tagint>(buf[m++]);
wfd_list[nlocal][n] = static_cast<float>(buf[m++]);
wf_list[nlocal][n] = static_cast<float>(buf[m++]);
degradation_ij[nlocal][n] = static_cast<float>(buf[m++]);
energy_per_bond[nlocal][n] = static_cast<float>(buf[m++]);
}
return m;
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_restart(int i, double *buf) {
int m = 0;
buf[m++] = 4 * npartner[i] + 2;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::unpack_restart(int nlocal, int nth) {
// ipage = NULL if being called from granular pair style init()
// skip to Nth set of extra values
// double **extra = atom->extra;
//
// int m = 0;
// for (int i = 0; i < nth; i++)
// m += static_cast<int>(extra[nlocal][m]);
// m++;
//
// // allocate new chunks from ipage,dpage for incoming values
//
// npartner[nlocal] = static_cast<int>(extra[nlocal][m++]);
// for (int n = 0; n < npartner[nlocal]; n++) {
// partner[nlocal][n] = static_cast<tagint>(extra[nlocal][m++]);
// }
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::maxsize_restart() {
// maxtouch_all = max # of touching partners across all procs
int maxtouch_all;
MPI_Allreduce(&maxpartner, &maxtouch_all, 1, MPI_INT, MPI_MAX, world);
return 4 * maxtouch_all + 2;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::size_restart(int nlocal) {
return 4 * npartner[nlocal] + 2;
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
int i, j, m;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
//printf("FixSMD_TLSPH_ReferenceConfiguration:::pack_forward_comm\n");
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
buf[m++] = vfrac[j];
buf[m++] = radius[j];
buf[m++] = defgrad0[i][0];
buf[m++] = defgrad0[i][1];
buf[m++] = defgrad0[i][2];
buf[m++] = defgrad0[i][3];
buf[m++] = defgrad0[i][4];
buf[m++] = defgrad0[i][5];
buf[m++] = defgrad0[i][6];
buf[m++] = defgrad0[i][7];
buf[m++] = defgrad0[i][8];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::unpack_forward_comm(int n, int first, double *buf) {
int i, m, last;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
vfrac[i] = buf[m++];
radius[i] = buf[m++];
defgrad0[i][0] = buf[m++];
defgrad0[i][1] = buf[m++];
defgrad0[i][2] = buf[m++];
defgrad0[i][3] = buf[m++];
defgrad0[i][4] = buf[m++];
defgrad0[i][5] = buf[m++];
defgrad0[i][6] = buf[m++];
defgrad0[i][7] = buf[m++];
defgrad0[i][8] = buf[m++];
}
}
/* ----------------------------------------------------------------------
routine for excluding bonds across a hardcoded slit crack
Note that everything is scaled by lattice constant l0 to avoid
numerical inaccuracies.
------------------------------------------------------------------------- */
bool FixSMD_TLSPH_ReferenceConfiguration::crack_exclude(int i, int j) {
double **x = atom->x;
double l0 = domain->lattice->xlattice;
// line between pair of atoms i,j
double x1 = x[i][0] / l0;
double y1 = x[i][1] / l0;
double x2 = x[j][0] / l0;
double y2 = x[j][1] / l0;
// hardcoded crack line
double x3 = -0.1 / l0;
double y3 = ((int) 1.0 / l0) + 0.5;
//printf("y3 = %f\n", y3);
double x4 = 0.1 / l0 - 1.0 + 0.1;
double y4 = y3;
bool retVal = DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4);
return !retVal;
//return 1;
}

View File

@ -0,0 +1,86 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* This file is based on the FixShearHistory class.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(SMD_TLSPH_NEIGHBORS,FixSMD_TLSPH_ReferenceConfiguration)
#else
#ifndef LMP_FIX_SMD_TLSPH_REFERENCE_H
#define LMP_FIX_SMD_TLSPH_REFERENCE_H
#include "fix.h"
#include "my_page.h"
namespace LAMMPS_NS {
class FixSMD_TLSPH_ReferenceConfiguration: public Fix {
friend class Neighbor;
friend class PairTlsph;
public:
FixSMD_TLSPH_ReferenceConfiguration(class LAMMPS *, int, char **);
~FixSMD_TLSPH_ReferenceConfiguration();
int setmask();
void init();
void setup(int);
void pre_exchange();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
int pack_restart(int, double *);
void unpack_restart(int, int);
int size_restart(int);
int maxsize_restart();
bool crack_exclude(int i, int j);
bool get_line_intersection(int i, int j);
protected:
int updateFlag; // flag to update reference configuration
int nmax;
int maxpartner;
int *npartner; // # of touching partners of each atom
tagint **partner; // global atom IDs for the partners
float **wfd_list, **wf_list, **energy_per_bond;
float **degradation_ij; // per-pair interaction degradation status
class Pair *pair;
};
}
#endif
#endif

View File

@ -0,0 +1,507 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Mike Parks (SNL), Ezwanur Rahman, J.T. Foster (UTSA)
------------------------------------------------------------------------- */
#include "math.h"
#include "fix_smd_wall_surface.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "comm.h"
#include "update.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "pair.h"
#include "lattice.h"
#include "memory.h"
#include "error.h"
#include <Eigen/Eigen>
#include <stdio.h>
#include <iostream>
#include "atom_vec.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace Eigen;
using namespace std;
#define DELTA 16384
#define EPSILON 1.0e-6
enum {
LAYOUT_UNIFORM, LAYOUT_NONUNIFORM, LAYOUT_TILED
};
// several files
/* ---------------------------------------------------------------------- */
FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
restart_global = 0;
restart_peratom = 0;
first = 1;
//atom->add_callback(0);
//atom->add_callback(1);
if (narg != 6)
error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface");
filename.assign(arg[3]);
wall_particle_type = force->inumeric(FLERR, arg[4]);
wall_molecule_id = force->inumeric(FLERR, arg[5]);
if (wall_molecule_id < 65535) {
error->one(FLERR, "wall molcule id must be >= 65535\n");
}
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename.c_str());
printf("fix smd/wall_surface has particle type %d \n", wall_particle_type);
printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
/* ---------------------------------------------------------------------- */
FixSMDWallSurface::~FixSMDWallSurface() {
// unregister this fix so atom class doesn't invoke it any more
//atom->delete_callback(id, 0);
//atom->delete_callback(id, 1);
}
/* ---------------------------------------------------------------------- */
int FixSMDWallSurface::setmask() {
int mask = 0;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDWallSurface::init() {
if (!first)
return;
}
/* ----------------------------------------------------------------------
For minimization: setup as with dynamics
------------------------------------------------------------------------- */
void FixSMDWallSurface::min_setup(int vflag) {
setup(vflag);
}
/* ----------------------------------------------------------------------
create initial list of neighbor partners via call to neighbor->build()
must be done in setup (not init) since fix init comes before neigh init
------------------------------------------------------------------------- */
void FixSMDWallSurface::setup(int vflag) {
if (!first)
return;
first = 0;
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// insures all data atoms will be owned even with round-off
int triclinic = domain->triclinic;
double epsilon[3];
if (triclinic)
epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
else {
epsilon[0] = domain->prd[0] * EPSILON;
epsilon[1] = domain->prd[1] * EPSILON;
epsilon[2] = domain->prd[2] * EPSILON;
}
if (triclinic == 0) {
sublo[0] = domain->sublo[0];
subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1];
subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2];
subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0];
subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1];
subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2];
subhi[2] = domain->subhi_lamda[2];
}
if (comm->layout != LAYOUT_TILED) {
if (domain->xperiodic) {
if (comm->myloc[0] == 0)
sublo[0] -= epsilon[0];
if (comm->myloc[0] == comm->procgrid[0] - 1)
subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->myloc[1] == 0)
sublo[1] -= epsilon[1];
if (comm->myloc[1] == comm->procgrid[1] - 1)
subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->myloc[2] == 0)
sublo[2] -= epsilon[2];
if (comm->myloc[2] == comm->procgrid[2] - 1)
subhi[2] += epsilon[2];
}
} else {
if (domain->xperiodic) {
if (comm->mysplit[0][0] == 0.0)
sublo[0] -= epsilon[0];
if (comm->mysplit[0][1] == 1.0)
subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->mysplit[1][0] == 0.0)
sublo[1] -= epsilon[1];
if (comm->mysplit[1][1] == 1.0)
subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->mysplit[2][0] == 0.0)
sublo[2] -= epsilon[2];
if (comm->mysplit[2][1] == 1.0)
subhi[2] += epsilon[2];
}
}
read_triangles(0);
}
/* ----------------------------------------------------------------------
function to determine number of values in a text line
------------------------------------------------------------------------- */
int FixSMDWallSurface::count_words(const char *line) {
int n = strlen(line) + 1;
char *copy;
memory->create(copy, n, "atom:copy");
strcpy(copy, line);
char *ptr;
if ((ptr = strchr(copy, '#')))
*ptr = '\0';
if (strtok(copy, " \t\n\r\f") == NULL) {
memory->destroy(copy);
return 0;
}
n = 1;
while (strtok(NULL, " \t\n\r\f"))
n++;
memory->destroy(copy);
return n;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
void FixSMDWallSurface::read_triangles(int pass) {
double coord[3];
int nlocal_previous = atom->nlocal;
int ilocal = nlocal_previous;
int m;
int me;
bigint natoms_previous = atom->natoms;
Vector3d *vert;
vert = new Vector3d[3];
Vector3d normal, center;
FILE *fp = fopen(filename.c_str(), "r");
if (fp == NULL) {
char str[128];
sprintf(str, "Cannot open file %s", filename.c_str());
error->one(FLERR, str);
}
MPI_Comm_rank(world, &me);
if (me == 0) {
if (screen) {
if (pass == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, " scanning triangle pairs ...\n");
} else {
fprintf(screen, " reading triangle pairs ...\n");
}
}
if (logfile) {
if (pass == 0) {
fprintf(logfile, " scanning triangle pairs ...\n");
} else {
fprintf(logfile, " reading triangle pairs ...\n");
}
}
}
char str[128];
char line[256];
char *retpointer;
char **values;
int nwords;
// read STL solid name
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading number of triangle pairs");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords < 1) {
sprintf(str, "first line of file is incorrect");
error->one(FLERR, str);
}
// values = new char*[nwords];
// values[0] = strtok(line, " \t\n\r\f");
// if (values[0] == NULL)
// error->all(FLERR, "Incorrect atom format in data file");
// for (m = 1; m < nwords; m++) {
// values[m] = strtok(NULL, " \t\n\r\f");
// if (values[m] == NULL)
// error->all(FLERR, "Incorrect atom format in data file");
// }
// delete[] values;
//
// if (comm->me == 0) {
// cout << "STL file contains solid body with name: " << values[1] << endl;
// }
// iterate over STL facets util end of body is reached
while (fgets(line, sizeof(line), fp)) { // read a line, should be the facet line
// evaluate facet line
nwords = count_words(line);
if (nwords != 5) {
//sprintf(str, "found end solid line");
//error->message(FLERR, str);
break;
} else {
// should be facet line
}
values = new char*[nwords];
values[0] = strtok(line, " \t\n\r\f");
if (values[0] == NULL)
error->all(FLERR, "Incorrect atom format in data file");
for (m = 1; m < nwords; m++) {
values[m] = strtok(NULL, " \t\n\r\f");
if (values[m] == NULL)
error->all(FLERR, "Incorrect atom format in data file");
}
normal << force->numeric(FLERR, values[2]), force->numeric(FLERR, values[3]), force->numeric(FLERR, values[4]);
//cout << "normal is " << normal << endl;
delete[] values;
// read outer loop line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading outer loop");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 2) {
sprintf(str, "error reading outer loop");
error->one(FLERR, str);
}
// read vertex lines
for (int k = 0; k < 3; k++) {
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading vertex line");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 4) {
sprintf(str, "error reading vertex line");
error->one(FLERR, str);
}
values = new char*[nwords];
values[0] = strtok(line, " \t\n\r\f");
if (values[0] == NULL)
error->all(FLERR, "Incorrect vertex line");
for (m = 1; m < nwords; m++) {
values[m] = strtok(NULL, " \t\n\r\f");
if (values[m] == NULL)
error->all(FLERR, "Incorrect vertex line");
}
vert[k] << force->numeric(FLERR, values[1]), force->numeric(FLERR, values[2]), force->numeric(FLERR, values[3]);
//cout << "vertex is " << vert[k] << endl;
//printf("%s %s %s\n", values[1], values[2], values[3]);
delete[] values;
//exit(1);
}
// read end loop line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading endloop");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 1) {
sprintf(str, "error reading endloop");
error->one(FLERR, str);
}
// read end facet line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading endfacet");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 1) {
sprintf(str, "error reading endfacet");
error->one(FLERR, str);
}
// now we have a normal and three vertices ... proceed with adding triangle
center = (vert[0] + vert[1] + vert[2]) / 3.0;
// cout << "center is " << center << endl;
double r1 = (center - vert[0]).norm();
double r2 = (center - vert[1]).norm();
double r3 = (center - vert[2]).norm();
double r = MAX(r1, r2);
r = MAX(r, r3);
/*
* if atom/molecule is in my subbox, create it
* ... use x0 to hold triangle normal.
* ... use smd_data_9 to hold the three vertices
* ... use x to hold triangle center
* ... radius is the mmaximal distance from triangle center to all vertices
*/
// printf("coord: %f %f %f\n", coord[0], coord[1], coord[2]);
// printf("sublo: %f %f %f\n", sublo[0], sublo[1], sublo[2]);
// printf("subhi: %f %f %f\n", subhi[0], subhi[1], subhi[2]);
//printf("ilocal = %d\n", ilocal);
if (center(0) >= sublo[0] && center(0) < subhi[0] && center(1) >= sublo[1] && center(1) < subhi[1] && center(2) >= sublo[2]
&& center(2) < subhi[2]) {
//printf("******* KERATIN nlocal=%d ***\n", nlocal);
coord[0] = center(0);
coord[1] = center(1);
coord[2] = center(2);
atom->avec->create_atom(wall_particle_type, coord);
/*
* need to initialize pointers to atom vec arrays here, because they could have changed
* due to calling grow() in create_atoms() above;
*/
int *mol = atom->molecule;
int *type = atom->type;
double *radius = atom->radius;
double *contact_radius = atom->contact_radius;
double **smd_data_9 = atom->smd_data_9;
double **x0 = atom->x0;
radius[ilocal] = r; //ilocal;
contact_radius[ilocal] = r; //ilocal;
mol[ilocal] = wall_molecule_id;
type[ilocal] = wall_particle_type;
x0[ilocal][0] = normal(0);
x0[ilocal][1] = normal(1);
x0[ilocal][2] = normal(2);
smd_data_9[ilocal][0] = vert[0](0);
smd_data_9[ilocal][1] = vert[0](1);
smd_data_9[ilocal][2] = vert[0](2);
smd_data_9[ilocal][3] = vert[1](0);
smd_data_9[ilocal][4] = vert[1](1);
smd_data_9[ilocal][5] = vert[1](2);
smd_data_9[ilocal][6] = vert[2](0);
smd_data_9[ilocal][7] = vert[2](1);
smd_data_9[ilocal][8] = vert[2](2);
ilocal++;
}
}
// set new total # of atoms and error check
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal, &atom->natoms, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (atom->natoms < 0 || atom->natoms >= MAXBIGINT)
error->all(FLERR, "Too many total atoms");
// add IDs for newly created atoms
// check that atom IDs are valid
if (atom->tag_enable)
atom->tag_extend();
atom->tag_check();
// create global mapping of atoms
// zero nghost in case are adding new atoms to existing atoms
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// print status
if (comm->me == 0) {
if (screen) {
printf("... fix smd/wall_surface finished reading triangulated surface\n");
fprintf(screen, "fix smd/wall_surface created " BIGINT_FORMAT " atoms\n", atom->natoms - natoms_previous);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
if (logfile) {
fprintf(logfile, "... fix smd/wall_surface finished reading triangulated surface\n");
fprintf(logfile, "fix smd/wall_surface created " BIGINT_FORMAT " atoms\n", atom->natoms - natoms_previous);
fprintf(logfile, ">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
delete[] vert;
fclose(fp);
}

View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(smd/wall_surface,FixSMDWallSurface)
#else
#ifndef LMP_FIX_SMD_WALL_SURFACE_H
#define LMP_FIX_SMD_WALL_SURFACE_H
#include "fix.h"
#include <iostream>
using namespace std;
namespace LAMMPS_NS {
class FixSMDWallSurface: public Fix {
public:
FixSMDWallSurface(class LAMMPS *, int, char **);
virtual ~FixSMDWallSurface();
int setmask();
void init();
void setup(int);
void min_setup(int);
int count_words(const char *line);
void read_triangles(int pass);
private:
int first; // flag for first time initialization
double sublo[3], subhi[3]; // epsilon-extended proc sub-box for adding atoms;
std::string filename;
int wall_particle_type;
int wall_molecule_id;
};
}
#endif
#endif

View File

@ -0,0 +1,385 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mike Parks (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "float.h"
#include "stdlib.h"
#include "string.h"
#include "pair_smd_hertz.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define SQRT2 1.414213562e0
/* ---------------------------------------------------------------------- */
PairHertz::PairHertz(LAMMPS *lmp) :
Pair(lmp) {
onerad_dynamic = onerad_frozen = maxrad_dynamic = maxrad_frozen = NULL;
bulkmodulus = NULL;
kn = NULL;
scale = 1.0;
}
/* ---------------------------------------------------------------------- */
PairHertz::~PairHertz() {
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(bulkmodulus);
memory->destroy(kn);
delete[] onerad_dynamic;
delete[] onerad_frozen;
delete[] maxrad_dynamic;
delete[] maxrad_frozen;
}
}
/* ---------------------------------------------------------------------- */
void PairHertz::compute(int eflag, int vflag) {
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz;
double rsq, r, evdwl, fpair;
int *ilist, *jlist, *numneigh, **firstneigh;
double rcut, r_geom, delta, ri, rj, dt_crit;
double *rmass = atom->rmass;
evdwl = 0.0;
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
double **f = atom->f;
double **x = atom->x;
double **x0 = atom->x0;
int *type = atom->type;
int nlocal = atom->nlocal;
double *radius = atom->contact_radius;
double *sph_radius = atom->radius;
double rcutSq;
double delx0, dely0, delz0, rSq0, sphCut;
int newton_pair = force->newton_pair;
int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
stable_time_increment = 1.0e22;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
ri = scale * radius[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
rj = scale * radius[j];
rcut = ri + rj;
rcutSq = rcut * rcut;
if (rsq < rcutSq) {
/*
* self contact option:
* if pair of particles was initially close enough to interact via a bulk continuum mechanism (e.g. SPH), exclude pair from contact forces.
* this approach should work well if no updates of the reference configuration are performed.
*/
if (itype == jtype) {
delx0 = x0[j][0] - x0[i][0];
dely0 = x0[j][1] - x0[i][1];
delz0 = x0[j][2] - x0[i][2];
if (periodic) {
domain->minimum_image(delx0, dely0, delz0);
}
rSq0 = delx0 * delx0 + dely0 * dely0 + delz0 * delz0; // initial distance
sphCut = sph_radius[i] + sph_radius[j];
if (rSq0 < sphCut * sphCut) {
rcut = 0.5 * rcut;
rcutSq = rcut * rcut;
if (rsq > rcutSq) {
continue;
}
}
}
r = sqrt(rsq);
//printf("hertz interaction, r=%f, cut=%f, h=%f\n", r, rcut, sqrt(rSq0));
// Hertzian short-range forces
delta = rcut - r; // overlap distance
r_geom = ri * rj / rcut;
//assuming poisson ratio = 1/4 for 3d
fpair = 1.066666667e0 * bulkmodulus[itype][jtype] * delta * sqrt(delta * r_geom); // units: N
evdwl = fpair * 0.4e0 * delta; // GCG 25 April: this expression conserves total energy
dt_crit = 3.14 * sqrt(0.5 * (rmass[i] + rmass[j]) / (fpair / delta));
stable_time_increment = MIN(stable_time_increment, dt_crit);
if (r > 2.0e-16) {
fpair /= r; // divide by r and multiply with non-normalized distance vector
} else {
fpair = 0.0;
}
/*
* contact viscosity -- needs to be done, see GRANULAR package for normal & shear damping
* for now: no damping and thus no viscous energy deltaE
*/
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
}
}
// double stable_time_increment_all = 0.0;
// MPI_Allreduce(&stable_time_increment, &stable_time_increment_all, 1, MPI_DOUBLE, MPI_MIN, world);
// if (comm->me == 0) {
// printf("stable time step for pair smd/hertz is %f\n", stable_time_increment_all);
// }
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairHertz::allocate() {
allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(bulkmodulus, n + 1, n + 1, "pair:kspring");
memory->create(kn, n + 1, n + 1, "pair:kn");
memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); // always needs to be allocated, even with granular neighborlist
onerad_dynamic = new double[n + 1];
onerad_frozen = new double[n + 1];
maxrad_dynamic = new double[n + 1];
maxrad_frozen = new double[n + 1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairHertz::settings(int narg, char **arg) {
if (narg != 1)
error->all(FLERR, "Illegal number of args for pair_style hertz");
scale = force->numeric(FLERR, arg[0]);
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("SMD/HERTZ CONTACT SETTINGS:\n");
printf("... effective contact radius is scaled by %f\n", scale);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairHertz::coeff(int narg, char **arg) {
if (narg != 3)
error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated)
allocate();
int ilo, ihi, jlo, jhi;
force->bounds(arg[0], atom->ntypes, ilo, ihi);
force->bounds(arg[1], atom->ntypes, jlo, jhi);
double bulkmodulus_one = atof(arg[2]);
// set short-range force constant
double kn_one = 0.0;
if (domain->dimension == 3) {
kn_one = (16. / 15.) * bulkmodulus_one; //assuming poisson ratio = 1/4 for 3d
} else {
kn_one = 0.251856195 * (2. / 3.) * bulkmodulus_one; //assuming poisson ratio = 1/3 for 2d
}
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
bulkmodulus[i][j] = bulkmodulus_one;
kn[i][j] = kn_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0)
error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairHertz::init_one(int i, int j) {
if (!allocated)
allocate();
if (setflag[i][j] == 0)
error->all(FLERR, "All pair coeffs are not set");
bulkmodulus[j][i] = bulkmodulus[i][j];
kn[j][i] = kn[i][j];
// cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
double cutoff = maxrad_dynamic[i] + maxrad_dynamic[j];
cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j]);
cutoff = MAX(cutoff, maxrad_dynamic[i] + maxrad_frozen[j]);
if (comm->me == 0) {
printf("cutoff for pair smd/hertz = %f\n", cutoff);
}
return cutoff;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairHertz::init_style() {
int i;
// error checks
if (!atom->contact_radius_flag)
error->all(FLERR, "Pair style smd/hertz requires atom style with contact_radius");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->gran = 1;
// set maxrad_dynamic and maxrad_frozen for each type
// include future Fix pour particles as dynamic
for (i = 1; i <= atom->ntypes; i++)
onerad_dynamic[i] = onerad_frozen[i] = 0.0;
double *radius = atom->radius;
int *type = atom->type;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]], radius[i]);
}
MPI_Allreduce(&onerad_dynamic[1], &maxrad_dynamic[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
MPI_Allreduce(&onerad_frozen[1], &maxrad_frozen[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
optional granular history list
------------------------------------------------------------------------- */
void PairHertz::init_list(int id, NeighList *ptr) {
if (id == 0)
list = ptr;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairHertz::memory_usage() {
return 0.0;
}
void *PairHertz::extract(const char *str, int &i) {
//printf("in PairTriSurf::extract\n");
if (strcmp(str, "smd/hertz/stable_time_increment_ptr") == 0) {
return (void *) &stable_time_increment;
}
return NULL;
}

View File

@ -0,0 +1,69 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(smd/hertz,PairHertz)
#else
#ifndef LMP_SMD_HERTZ_H
#define LMP_SMD_HERTZ_H
#include "pair.h"
namespace LAMMPS_NS {
class PairHertz : public Pair {
public:
PairHertz(class LAMMPS *);
virtual ~PairHertz();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void init_list(int, class NeighList *);
virtual double memory_usage();
void *extract(const char *, int &);
protected:
double **bulkmodulus;
double **kn;
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;
double scale;
double stable_time_increment; // stable time step size
void allocate();
};
}
#endif
#endif

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,229 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(smd/tlsph,PairTlsph)
#else
#ifndef LMP_TLSPH_NEW_H
#define LMP_TLSPH_NEW_H
#include "pair.h"
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <map>
using namespace std;
using namespace Eigen;
namespace LAMMPS_NS {
class PairTlsph: public Pair {
public:
PairTlsph(class LAMMPS *);
virtual ~PairTlsph();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void init_list(int, class NeighList *);
void write_restart_settings(FILE *) {
}
void read_restart_settings(FILE *) {
}
virtual double memory_usage();
void compute_shape_matrix(void);
void material_model(void);
void *extract(const char *, int &);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
void AssembleStress();
void PreCompute();
void ComputeForces(int eflag, int vflag);
void effective_longitudinal_modulus(const int itype, const double dt, const double d_iso, const double p_rate,
const Matrix3d d_dev, const Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff);
void ComputePressure(const int i, const double rho, const double mass_specific_energy, const double vol_specific_energy,
const double pInitial, const double d_iso, double &pFinal, double &p_rate);
void ComputeStressDeviator(const int i, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev,
Matrix3d &sigma_dev_rate, double &plastic_strain_increment);
void ComputeDamage(const int i, const Matrix3d strain, const Matrix3d sigmaFinal, Matrix3d &sigma_damaged);
protected:
void allocate();
char *suffix;
/*
* per-type arrays
*/
int *strengthModel, *eos;
double *onerad_dynamic, *onerad_frozen, *maxrad_dynamic, *maxrad_frozen;
/*
* per atom arrays
*/
Matrix3d *K, *PK1, *Fdot, *Fincr;
Matrix3d *R; // rotation matrix
Matrix3d *FincrInv;
Matrix3d *D, *W; // strain rate and spin tensor
Vector3d *smoothVelDifference;
Matrix3d *CauchyStress;
double *detF, *particle_dt;
double *hourglass_error;
int *numNeighsRefConfig;
int nmax; // max number of atoms on this proc
double hMin; // minimum kernel radius for two particles
double dtCFL;
double dtRelative; // relative velocity of two particles, divided by sound speed
int updateFlag;
double update_threshold; // updateFlage is set to one if the relative displacement of a pair exceeds update_threshold
double cut_comm;
enum {
UPDATE_NONE = 5000, UPDATE_CONSTANT_THRESHOLD = 5001, UPDATE_PAIRWISE_RATIO = 5002,
};
enum {
LINEAR_DEFGRAD = 0,
STRENGTH_LINEAR = 1,
STRENGTH_LINEAR_PLASTIC = 2,
STRENGTH_JOHNSON_COOK = 3,
STRENGTH_NONE = 4,
EOS_LINEAR = 5,
EOS_SHOCK = 6,
EOS_POLYNOMIAL = 7,
EOS_NONE = 8,
REFERENCE_DENSITY = 9,
YOUNGS_MODULUS = 10,
POISSON_RATIO = 11,
HOURGLASS_CONTROL_AMPLITUDE = 12,
HEAT_CAPACITY = 13,
LAME_LAMBDA = 14,
SHEAR_MODULUS = 15,
M_MODULUS = 16,
SIGNAL_VELOCITY = 17,
BULK_MODULUS = 18,
VISCOSITY_Q1 = 19,
VISCOSITY_Q2 = 20,
YIELD_STRESS = 21,
FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD = 22,
JC_A = 23,
JC_B = 24,
JC_a = 25,
JC_C = 26,
JC_epdot0 = 27,
JC_T0 = 28,
JC_Tmelt = 29,
JC_M = 30,
EOS_SHOCK_C0 = 31,
EOS_SHOCK_S = 32,
EOS_SHOCK_GAMMA = 33,
HARDENING_PARAMETER = 34,
FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD = 35,
FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD = 36,
FAILURE_MAX_PAIRWISE_STRAIN_THRESHOLD = 37,
EOS_POLYNOMIAL_C0 = 38,
EOS_POLYNOMIAL_C1 = 39,
EOS_POLYNOMIAL_C2 = 40,
EOS_POLYNOMIAL_C3 = 41,
EOS_POLYNOMIAL_C4 = 42,
EOS_POLYNOMIAL_C5 = 43,
EOS_POLYNOMIAL_C6 = 44,
FAILURE_JC_D1 = 45,
FAILURE_JC_D2 = 46,
FAILURE_JC_D3 = 47,
FAILURE_JC_D4 = 48,
FAILURE_JC_EPDOT0 = 49,
CRITICAL_ENERGY_RELEASE_RATE = 50,
MAX_KEY_VALUE = 51
};
struct failure_types { // this is defined per type and determines which failure/damage model is active
bool failure_none;
bool failure_max_principal_strain;
bool failure_max_principal_stress;
bool failure_max_plastic_strain;
bool failure_johnson_cook;
bool failure_max_pairwise_strain;
bool integration_point_wise; // true if failure model applies to stress/strain state of integration point
bool failure_energy_release_rate;
failure_types() {
failure_none = true;
failure_max_principal_strain = false;
failure_max_principal_stress = false;
failure_max_plastic_strain = false;
failure_johnson_cook = false;
failure_max_pairwise_strain = false;
integration_point_wise = false;
failure_energy_release_rate = false;
//printf("constructed failure type\n");
}
};
failure_types *failureModel;
int ifix_tlsph;
int update_method;
class FixSMD_TLSPH_ReferenceConfiguration *fix_tlsph_reference_configuration;
private:
double **Lookup; // holds per-type material parameters for the quantities defined in enum statement above.
bool first; // if first is true, do not perform any computations, beacuse reference configuration is not ready yet.
};
}
#endif
#endif
/*
* materialCoeffs array for EOS parameters:
* 1: rho0
*
*
* materialCoeffs array for strength parameters:
*
* Common
* 10: maximum strain threshold for damage model
* 11: maximum stress threshold for damage model
*
* Linear Plasticity model:
* 12: plastic yield stress
*
*
* Blei: rho = 11.34e-6, c0=2000, s=1.46, Gamma=2.77
* Stahl 1403: rho = 7.86e-3, c=4569, s=1.49, Gamma=2.17
*/

View File

@ -0,0 +1,846 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mike Parks (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "float.h"
#include "stdlib.h"
#include "string.h"
#include "pair_smd_triangulated_surface.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include <Eigen/Eigen>
#include <stdio.h>
#include <iostream>
using namespace std;
using namespace LAMMPS_NS;
using namespace Eigen;
#define SQRT2 1.414213562e0
/* ---------------------------------------------------------------------- */
PairTriSurf::PairTriSurf(LAMMPS *lmp) :
Pair(lmp) {
onerad_dynamic = onerad_frozen = maxrad_dynamic = maxrad_frozen = NULL;
bulkmodulus = NULL;
kn = NULL;
scale = 1.0;
}
/* ---------------------------------------------------------------------- */
PairTriSurf::~PairTriSurf() {
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(bulkmodulus);
memory->destroy(kn);
delete[] onerad_dynamic;
delete[] onerad_frozen;
delete[] maxrad_dynamic;
delete[] maxrad_frozen;
}
}
/* ---------------------------------------------------------------------- */
void PairTriSurf::compute(int eflag, int vflag) {
int i, j, ii, jj, inum, jnum, itype, jtype;
double rsq, r, evdwl, fpair;
int *ilist, *jlist, *numneigh, **firstneigh;
double rcut, r_geom, delta, r_tri, r_particle, touch_distance, dt_crit;
int tri, particle;
Vector3d normal, x1, x2, x3, x4, x13, x23, x43, w, cp, x4cp, vnew, v_old;
;
Vector3d xi, x_center, dx;
Matrix2d C;
Vector2d w2d, rhs;
evdwl = 0.0;
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
int *mol = atom->molecule;
double **f = atom->f;
double **smd_data_9 = atom->smd_data_9;
double **x = atom->x;
double **x0 = atom->x0;
double **v = atom->v;
double *rmass = atom->rmass;
int *type = atom->type;
int nlocal = atom->nlocal;
double *radius = atom->contact_radius;
double rcutSq;
Vector3d offset;
int newton_pair = force->newton_pair;
int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
int max_neighs = 0;
stable_time_increment = 1.0e22;
// loop over neighbors of my atoms using a half neighbor list
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
max_neighs = MAX(max_neighs, jnum);
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
/*
* decide which one of i, j is triangle and which is particle
*/
if ((mol[i] < 65535) && (mol[j] >= 65535)) {
particle = i;
tri = j;
} else if ((mol[j] < 65535) && (mol[i] >= 65535)) {
particle = j;
tri = i;
} else {
error->one(FLERR, "unknown case");
}
//x_center << x[tri][0], x[tri][1], x[tri][2]; // center of triangle
x_center(0) = x[tri][0];
x_center(1) = x[tri][1];
x_center(2) = x[tri][2];
//x4 << x[particle][0], x[particle][1], x[particle][2];
x4(0) = x[particle][0];
x4(1) = x[particle][1];
x4(2) = x[particle][2];
dx = x_center - x4; //
if (periodic) {
domain->minimum_image(dx(0), dx(1), dx(2));
}
rsq = dx.squaredNorm();
r_tri = scale * radius[tri];
r_particle = scale * radius[particle];
rcut = r_tri + r_particle;
rcutSq = rcut * rcut;
//printf("type i=%d, type j=%d, r=%f, ri=%f, rj=%f\n", itype, jtype, sqrt(rsq), ri, rj);
if (rsq < rcutSq) {
/*
* gather triangle information
*/
normal(0) = x0[tri][0];
normal(1) = x0[tri][1];
normal(2) = x0[tri][2];
/*
* distance check: is particle closer than its radius to the triangle plane?
*/
if (fabs(dx.dot(normal)) < radius[particle]) {
/*
* get other two triangle vertices
*/
x1(0) = smd_data_9[tri][0];
x1(1) = smd_data_9[tri][1];
x1(2) = smd_data_9[tri][2];
x2(0) = smd_data_9[tri][3];
x2(1) = smd_data_9[tri][4];
x2(2) = smd_data_9[tri][5];
x3(0) = smd_data_9[tri][6];
x3(1) = smd_data_9[tri][7];
x3(2) = smd_data_9[tri][8];
PointTriangleDistance(x4, x1, x2, x3, cp, r);
/*
* distance to closest point
*/
x4cp = x4 - cp;
/*
* flip normal to point in direction of x4cp
*/
if (x4cp.dot(normal) < 0.0) {
normal *= -1.0;
}
/*
* penalty force pushes particle away from triangle
*/
if (r < 1.0 * radius[particle]) {
delta = radius[particle] - r; // overlap distance
r_geom = radius[particle];
fpair = 1.066666667e0 * bulkmodulus[itype][jtype] * delta * sqrt(delta * r_geom);
dt_crit = 3.14 * sqrt(rmass[particle] / (fpair / delta));
stable_time_increment = MIN(stable_time_increment, dt_crit);
evdwl = r * fpair * 0.4e0 * delta; // GCG 25 April: this expression conserves total energy
fpair /= (r + 1.0e-2 * radius[particle]); // divide by r + softening and multiply with non-normalized distance vector
if (particle < nlocal) {
f[particle][0] += x4cp(0) * fpair;
f[particle][1] += x4cp(1) * fpair;
f[particle][2] += x4cp(2) * fpair;
}
if (tri < nlocal) {
f[tri][0] -= x4cp(0) * fpair;
f[tri][1] -= x4cp(1) * fpair;
f[tri][2] -= x4cp(2) * fpair;
}
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, x4cp(0), x4cp(1), x4cp(2));
}
}
/*
* if particle comes too close to triangle, reflect its velocity and explicitely move it away
*/
touch_distance = 1.0 * radius[particle];
if (r < touch_distance) {
/*
* reflect velocity if it points toward triangle
*/
normal = x4cp / r;
//v_old << v[particle][0], v[particle][1], v[particle][2];
v_old(0) = v[particle][0];
v_old(1) = v[particle][1];
v_old(2) = v[particle][2];
if (v_old.dot(normal) < 0.0) {
//printf("flipping velocity\n");
vnew = 1.0 * (-2.0 * v_old.dot(normal) * normal + v_old);
v[particle][0] = vnew(0);
v[particle][1] = vnew(1);
v[particle][2] = vnew(2);
}
//printf("moving particle on top of triangle\n");
x[particle][0] = cp(0) + touch_distance * normal(0);
x[particle][1] = cp(1) + touch_distance * normal(1);
x[particle][2] = cp(2) + touch_distance * normal(2);
}
}
}
}
}
// int max_neighs_all = 0;
// MPI_Allreduce(&max_neighs, &max_neighs_all, 1, MPI_INT, MPI_MAX, world);
// if (comm->me == 0) {
// printf("max. neighs in tri pair is %d\n", max_neighs_all);
// }
//
// double stable_time_increment_all = 0.0;
// MPI_Allreduce(&stable_time_increment, &stable_time_increment_all, 1, MPI_DOUBLE, MPI_MIN, world);
// if (comm->me == 0) {
// printf("stable time step tri pair is %f\n", stable_time_increment_all);
// }
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairTriSurf::allocate() {
allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(bulkmodulus, n + 1, n + 1, "pair:kspring");
memory->create(kn, n + 1, n + 1, "pair:kn");
memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); // always needs to be allocated, even with granular neighborlist
onerad_dynamic = new double[n + 1];
onerad_frozen = new double[n + 1];
maxrad_dynamic = new double[n + 1];
maxrad_frozen = new double[n + 1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairTriSurf::settings(int narg, char **arg) {
if (narg != 1)
error->all(FLERR, "Illegal number of args for pair_style smd/tri_surface");
scale = force->numeric(FLERR, arg[0]);
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("SMD/TRI_SURFACE CONTACT SETTINGS:\n");
printf("... effective contact radius is scaled by %f\n", scale);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairTriSurf::coeff(int narg, char **arg) {
if (narg != 3)
error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated)
allocate();
int ilo, ihi, jlo, jhi;
force->bounds(arg[0], atom->ntypes, ilo, ihi);
force->bounds(arg[1], atom->ntypes, jlo, jhi);
double bulkmodulus_one = atof(arg[2]);
// set short-range force constant
double kn_one = 0.0;
if (domain->dimension == 3) {
kn_one = (16. / 15.) * bulkmodulus_one; //assuming poisson ratio = 1/4 for 3d
} else {
kn_one = 0.251856195 * (2. / 3.) * bulkmodulus_one; //assuming poisson ratio = 1/3 for 2d
}
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
bulkmodulus[i][j] = bulkmodulus_one;
kn[i][j] = kn_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0)
error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairTriSurf::init_one(int i, int j) {
if (!allocated)
allocate();
if (setflag[i][j] == 0)
error->all(FLERR, "All pair coeffs are not set");
bulkmodulus[j][i] = bulkmodulus[i][j];
kn[j][i] = kn[i][j];
// cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
double cutoff = maxrad_dynamic[i] + maxrad_dynamic[j];
cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j]);
cutoff = MAX(cutoff, maxrad_dynamic[i] + maxrad_frozen[j]);
if (comm->me == 0) {
printf("cutoff for pair smd/smd/tri_surface = %f\n", cutoff);
}
return cutoff;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairTriSurf::init_style() {
int i;
// error checks
if (!atom->contact_radius_flag)
error->all(FLERR, "Pair style smd/smd/tri_surface requires atom style with contact_radius");
// old: half list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->gran = 1;
// need a full neighbor list
// int irequest = neighbor->request(this);
// neighbor->requests[irequest]->half = 0;
// neighbor->requests[irequest]->full = 1;
// set maxrad_dynamic and maxrad_frozen for each type
// include future Fix pour particles as dynamic
for (i = 1; i <= atom->ntypes; i++)
onerad_dynamic[i] = onerad_frozen[i] = 0.0;
double *radius = atom->radius;
int *type = atom->type;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]], radius[i]);
}
MPI_Allreduce(&onerad_dynamic[1], &maxrad_dynamic[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
MPI_Allreduce(&onerad_frozen[1], &maxrad_frozen[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
optional granular history list
------------------------------------------------------------------------- */
void PairTriSurf::init_list(int id, NeighList *ptr) {
if (id == 0)
list = ptr;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairTriSurf::memory_usage() {
return 0.0;
}
/*
* distance between triangle and point
*/
/*
function [dist,PP0] = pointTriangleDistance(TRI,P)
% calculate distance between a point and a triangle in 3D
% SYNTAX
% dist = pointTriangleDistance(TRI,P)
% [dist,PP0] = pointTriangleDistance(TRI,P)
%
% DESCRIPTION
% Calculate the distance of a given point P from a triangle TRI.
% Point P is a row vector of the form 1x3. The triangle is a matrix
% formed by three rows of points TRI = [P1;P2;P3] each of size 1x3.
% dist = pointTriangleDistance(TRI,P) returns the distance of the point P
% to the triangle TRI.
% [dist,PP0] = pointTriangleDistance(TRI,P) additionally returns the
% closest point PP0 to P on the triangle TRI.
%
% Author: Gwendolyn Fischer
% Release: 1.0
% Release date: 09/02/02
% Release: 1.1 Fixed Bug because of normalization
% Release: 1.2 Fixed Bug because of typo in region 5 20101013
% Release: 1.3 Fixed Bug because of typo in region 2 20101014
% Possible extention could be a version tailored not to return the distance
% and additionally the closest point, but instead return only the closest
% point. Could lead to a small speed gain.
% Example:
% %% The Problem
% P0 = [0.5 -0.3 0.5];
%
% P1 = [0 -1 0];
% P2 = [1 0 0];
% P3 = [0 0 0];
%
% vertices = [P1; P2; P3];
% faces = [1 2 3];
%
% %% The Engine
% [dist,PP0] = pointTriangleDistance([P1;P2;P3],P0);
%
% %% Visualization
% [x,y,z] = sphere(20);
% x = dist*x+P0(1);
% y = dist*y+P0(2);
% z = dist*z+P0(3);
%
% figure
% hold all
% patch('Vertices',vertices,'Faces',faces,'FaceColor','r','FaceAlpha',0.8);
% plot3(P0(1),P0(2),P0(3),'b*');
% plot3(PP0(1),PP0(2),PP0(3),'*g')
% surf(x,y,z,'FaceColor','b','FaceAlpha',0.3)
% view(3)
% The algorithm is based on
% "David Eberly, 'Distance Between Point and Triangle in 3D',
% Geometric Tools, LLC, (1999)"
% http:\\www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
%
% ^t
% \ |
% \reg2|
% \ |
% \ |
% \ |
% \|
% *P2
% |\
% | \
% reg3 | \ reg1
% | \
% |reg0\
% | \
% | \ P1
% -------*-------*------->s
% |P0 \
% reg4 | reg5 \ reg6
*/
//void PairTriSurf::PointTriangleDistance(const Vector3d P, const Vector3d TRI1, const Vector3d TRI2, const Vector3d TRI3,
// Vector3d &CP, double &dist) {
//
// Vector3d B, E0, E1, D;
// double a, b, c, d, e, f;
// double det, s, t, sqrDistance, tmp0, tmp1, numer, denom, invDet;
//
// // rewrite triangle in normal form
// B = TRI1;
// E0 = TRI2 - B;
// E1 = TRI3 - B;
//
// D = B - P;
// a = E0.dot(E0);
// b = E0.dot(E1);
// c = E1.dot(E1);
// d = E0.dot(D);
// e = E1.dot(D);
// f = D.dot(D);
//
// det = a * c - b * b;
// //% do we have to use abs here?
// s = b * e - c * d;
// t = b * d - a * e;
//
// //% Terible tree of conditionals to determine in which region of the diagram
// //% shown above the projection of the point into the triangle-plane lies.
// if ((s + t) <= det) {
// if (s < 0) {
// if (t < 0) {
// // %region4
// if (d < 0) {
// t = 0;
// if (-d >= a) {
// s = 1;
// sqrDistance = a + 2 * d + f;
// } else {
// s = -d / a;
// sqrDistance = d * s + f;
// }
// } else {
// s = 0;
// if (e >= 0) {
// t = 0;
// sqrDistance = f;
// } else {
// if (-e >= c) {
// t = 1;
// sqrDistance = c + 2 * e + f;
// } else {
// t = -e / c;
// sqrDistance = e * t + f;
// }
// }
// }
// // end % of region 4
// } else {
// // % region 3
// s = 0;
// if (e >= 0) {
// t = 0;
// sqrDistance = f;
// } else {
// if (-e >= c) {
// t = 1;
// sqrDistance = c + 2 * e + f;
// } else {
// t = -e / c;
// sqrDistance = e * t + f;
// }
// }
// }
// // end of region 3
// } else {
// if (t < 0) {
// //% region 5
// t = 0;
// if (d >= 0) {
// s = 0;
// sqrDistance = f;
// } else {
// if (-d >= a) {
// s = 1;
// sqrDistance = a + 2 * d + f;
// } else {
// s = -d / a;
// sqrDistance = d * s + f;
// }
// }
// } else {
// // region 0
// invDet = 1 / det;
// s = s * invDet;
// t = t * invDet;
// sqrDistance = s * (a * s + b * t + 2 * d) + t * (b * s + c * t + 2 * e) + f;
// }
// }
// } else {
// if (s < 0) {
// // % region 2
// tmp0 = b + d;
// tmp1 = c + e;
// if (tmp1 > tmp0) { //% minimum on edge s+t=1
// numer = tmp1 - tmp0;
// denom = a - 2 * b + c;
// if (numer >= denom) {
// s = 1;
// t = 0;
// sqrDistance = a + 2 * d + f;
// } else {
// s = numer / denom;
// t = 1 - s;
// sqrDistance = s * (a * s + b * t + 2 * d) + t * (b * s + c * t + 2 * e) + f;
// }
// } else
// // % minimum on edge s=0
// s = 0;
// if (tmp1 <= 0) {
// t = 1;
// sqrDistance = c + 2 * e + f;
// } else {
// if (e >= 0) {
// t = 0;
// sqrDistance = f;
// } else {
// t = -e / c;
// sqrDistance = e * t + f;
// }
// }
// } //end % of region 2
// else {
// if (t < 0) {
// // %region6
// tmp0 = b + e;
// tmp1 = a + d;
// if (tmp1 > tmp0) {
// numer = tmp1 - tmp0;
// denom = a - 2 * b + c;
// if (numer >= denom) {
// t = 1;
// s = 0;
// sqrDistance = c + 2 * e + f;
// } else {
// t = numer / denom;
// s = 1 - t;
// sqrDistance = s * (a * s + b * t + 2 * d) + t * (b * s + c * t + 2 * e) + f;
// }
// } else {
// t = 0;
// if (tmp1 <= 0) {
// s = 1;
// sqrDistance = a + 2 * d + f;
// } else {
// if (d >= 0) {
// s = 0;
// sqrDistance = f;
// } else {
// s = -d / a;
// sqrDistance = d * s + f;
// }
// }
// } // % end region 6
// } else {
// //% region 1
// numer = c + e - b - d;
// if (numer <= 0) {
// s = 0;
// t = 1;
// sqrDistance = c + 2 * e + f;
// } else {
// denom = a - 2 * b + c;
// if (numer >= denom) {
// s = 1;
// t = 0;
// sqrDistance = a + 2 * d + f;
// } else {
// s = numer / denom;
// t = 1 - s;
// sqrDistance = s * (a * s + b * t + 2 * d) + t * (b * s + c * t + 2 * e) + f;
// }
// } //% end of region 1
// }
// }
// }
//
// // % account for numerical round-off error
// if (sqrDistance < 0) {
// sqrDistance = 0;
// }
//
// dist = sqrt(sqrDistance);
//
// // closest point
// CP = B + s * E0 + t * E1;
//
//}
/*
* % The algorithm is based on
% "David Eberly, 'Distance Between Point and Triangle in 3D',
% Geometric Tools, LLC, (1999)"
% http:\\www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
*/
void PairTriSurf::PointTriangleDistance(const Vector3d sourcePosition, const Vector3d TRI0, const Vector3d TRI1,
const Vector3d TRI2, Vector3d &CP, double &dist) {
Vector3d edge0 = TRI1 - TRI0;
Vector3d edge1 = TRI2 - TRI0;
Vector3d v0 = TRI0 - sourcePosition;
double a = edge0.dot(edge0);
double b = edge0.dot(edge1);
double c = edge1.dot(edge1);
double d = edge0.dot(v0);
double e = edge1.dot(v0);
double det = a * c - b * b;
double s = b * e - c * d;
double t = b * d - a * e;
if (s + t < det) {
if (s < 0.f) {
if (t < 0.f) {
if (d < 0.f) {
s = clamp(-d / a, 0.f, 1.f);
t = 0.f;
} else {
s = 0.f;
t = clamp(-e / c, 0.f, 1.f);
}
} else {
s = 0.f;
t = clamp(-e / c, 0.f, 1.f);
}
} else if (t < 0.f) {
s = clamp(-d / a, 0.f, 1.f);
t = 0.f;
} else {
float invDet = 1.f / det;
s *= invDet;
t *= invDet;
}
} else {
if (s < 0.f) {
float tmp0 = b + d;
float tmp1 = c + e;
if (tmp1 > tmp0) {
float numer = tmp1 - tmp0;
float denom = a - 2 * b + c;
s = clamp(numer / denom, 0.f, 1.f);
t = 1 - s;
} else {
t = clamp(-e / c, 0.f, 1.f);
s = 0.f;
}
} else if (t < 0.f) {
if (a + d > b + e) {
float numer = c + e - b - d;
float denom = a - 2 * b + c;
s = clamp(numer / denom, 0.f, 1.f);
t = 1 - s;
} else {
s = clamp(-e / c, 0.f, 1.f);
t = 0.f;
}
} else {
float numer = c + e - b - d;
float denom = a - 2 * b + c;
s = clamp(numer / denom, 0.f, 1.f);
t = 1.f - s;
}
}
CP = TRI0 + s * edge0 + t * edge1;
dist = (CP - sourcePosition).norm();
}
double PairTriSurf::clamp(const double a, const double min, const double max) {
if (a < min) {
return min;
} else if (a > max) {
return max;
} else {
return a;
}
}
void *PairTriSurf::extract(const char *str, int &i) {
//printf("in PairTriSurf::extract\n");
if (strcmp(str, "smd/tri_surface/stable_time_increment_ptr") == 0) {
return (void *) &stable_time_increment;
}
return NULL;
}

View File

@ -0,0 +1,75 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(smd/tri_surface,PairTriSurf)
#else
#ifndef LMP_SMD_TRI_SURFACE_H
#define LMP_SMD_TRI_SURFACE_H
#include "pair.h"
#include <Eigen/Eigen>
using namespace Eigen;
namespace LAMMPS_NS {
class PairTriSurf : public Pair {
public:
PairTriSurf(class LAMMPS *);
virtual ~PairTriSurf();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void init_list(int, class NeighList *);
virtual double memory_usage();
void PointTriangleDistance(const Vector3d P, const Vector3d TRI1, const Vector3d TRI2, const Vector3d TRI3,
Vector3d &CP, double &dist);
double clamp(const double a, const double min, const double max);
void *extract(const char *, int &);
protected:
double **bulkmodulus;
double **kn;
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;
double scale;
double stable_time_increment; // stable time step size
void allocate();
};
}
#endif
#endif

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,140 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(smd/ulsph,PairULSPH)
#else
#ifndef LMP_ULSPH_H
#define LMP_ULSPH_H
#include "pair.h"
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/SVD>
using namespace Eigen;
using namespace std;
namespace LAMMPS_NS {
class PairULSPH: public Pair {
public:
PairULSPH(class LAMMPS *);
virtual ~PairULSPH();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void init_list(int, class NeighList *);
virtual double memory_usage();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
void AssembleStressTensor();
void *extract(const char *, int &);
void PreCompute();
void PreCompute_DensitySummation();
double effective_shear_modulus(const Matrix3d d_dev, const Matrix3d deltaStressDev, const double dt, const int itype);
Vector3d ComputeHourglassForce(const int i, const int itype, const int j, const int jtype, const Vector3d dv,
const Vector3d xij, const Vector3d g, const double c_ij, const double mu_ij, const double rho_ij);
protected:
double *c0_type; // reference speed of sound defined per particle type
double *rho0; // reference mass density per type
double *Q1; // linear artificial viscosity coeff
int *eos, *viscosity, *strength; // eos and strength material models
double **artificial_pressure; // true/false: use Monaghan's artificial pressure correction?
double **artificial_stress; // artificial stress amplitude
double *onerad_dynamic, *onerad_frozen;
double *maxrad_dynamic, *maxrad_frozen;
void allocate();
int nmax; // max number of atoms on this proc
int *numNeighs;
Matrix3d *K;
double *shepardWeight, *c0, *rho;
Vector3d *smoothVel;
Matrix3d *stressTensor, *L, *F;
double dtCFL;
private:
// enumerate EOSs. MUST BE IN THE RANGE [1000, 2000)
enum {
EOS_LINEAR = 1000, EOS_PERFECT_GAS = 1001, EOS_TAIT = 1002,
};
// enumerate physical viscosity models. MUST BE IN THE RANGE [2000, 3000)
enum {
VISCOSITY_NEWTON = 2000
};
// enumerate strength models. MUST BE IN THE RANGE [3000, 4000)
enum {
STRENGTH_LINEAR = 3000, STRENGTH_LINEAR_PLASTIC = 3001
};
// enumerate some quantitities and associate these with integer values such that they can be used for lookup in an array structure
enum {
NONE = 0,
BULK_MODULUS = 1,
HOURGLASS_CONTROL_AMPLITUDE = 2,
EOS_TAIT_EXPONENT = 3,
REFERENCE_SOUNDSPEED = 4,
REFERENCE_DENSITY = 5,
EOS_PERFECT_GAS_GAMMA = 6,
SHEAR_MODULUS = 7,
YIELD_STRENGTH = 8,
YOUNGS_MODULUS = 9,
POISSON_RATIO = 10,
LAME_LAMBDA = 11,
HEAT_CAPACITY = 12,
M_MODULUS = 13,
HARDENING_PARAMETER = 14,
VISCOSITY_MU = 15,
MAX_KEY_VALUE = 16
};
double **Lookup; // holds per-type material parameters for the quantities defined in enum statement above.
bool velocity_gradient_required;
int updateFlag; // indicates if any relative particle pair movement is significant compared to smoothing length
bool density_summation, density_continuity, velocity_gradient, gradient_correction_flag;
double *effm;
};
}
#endif
#endif

146
src/USER-SMD/smd_kernels.h Normal file
View File

@ -0,0 +1,146 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
#ifndef SMD_KERNEL_FUNCTIONS_H_
#define SMD_KERNEL_FUNCTIONS_H_
namespace SMD_Kernels {
static inline double Kernel_Wendland_Quintic_NotNormalized(const double r, const double h) {
if (r < h) {
double q = 2.0 * r / h;
return pow(1.0 - 0.5 * q, 4) * (2.0 * q + 1.0);
} else {
return 0.0;
}
}
static inline double Kernel_Cubic_Spline(const double r, const double h) {
double q = 2.0 * r / h;
if (q > 2.0) {
return 0.0;
} else if ((q <= 2.0) && (q > 1.0)) {
return pow(2.0 - q, 3.0) / 6.0;
} else if ((q >= 0.0) && (q <= 1.0)) {
return 2. / 3. - q * q + 0.5 * q * q * q;
} else {
return 0.0;
}
}
static inline double Kernel_Barbara(const double r, const double h) {
double arg = (1.570796327 * (r + h)) / h;
double hsq = h * h;
//wf = (1.680351548 * (cos(arg) + 1.)) / hsq;
return -2.639490040 * sin(arg) / (hsq * h);
}
static inline void spiky_kernel_and_derivative(const double h, const double r, const int dimension, double &wf, double &wfd) {
/*
* Spiky kernel
*/
if (r > h) {
printf("r=%f > h=%f in Spiky kernel\n", r, h);
wf = wfd = 0.0;
return;
}
double hr = h - r; // [m]
if (dimension == 2) {
double n = 0.3141592654e0 * h * h * h * h * h; // [m^5]
wfd = -3.0e0 * hr * hr / n; // [m*m/m^5] = [1/m^3] ==> correct for dW/dr in 2D
wf = -0.333333333333e0 * hr * wfd; // [m/m^3] ==> [1/m^2] correct for W in 2D
} else {
wfd = -14.0323944878e0 * hr * hr / (h * h * h * h * h * h); // [1/m^4] ==> correct for dW/dr in 3D
wf = -0.333333333333e0 * hr * wfd; // [m/m^4] ==> [1/m^3] correct for W in 3D
}
// alternative formulation
// double hr = h - r;
//
// /*
// * Spiky kernel
// */
//
// if (domain->dimension == 2) {
// double h5 = h * h * h * h * h;
// wf = 3.183098861e0 * hr * hr * hr / h5;
// wfd = -9.549296583 * hr * hr / h5;
//
// } else {
// double h6 = h * h * h * h * h * h;
// wf = 4.774648292 * hr * hr * hr / h6;
// wfd = -14.32394487 * hr * hr / h6;
// }
// }
}
static inline void barbara_kernel_and_derivative(const double h, const double r, const int dimension, double &wf, double &wfd) {
/*
* Barbara kernel
*/
double arg = (1.570796327 * (r + h)) / h;
double hsq = h * h;
if (r > h) {
printf("r = %f > h = %f in barbara kernel function\n", r, h);
exit(1);
//wf = wfd = 0.0;
//return;
}
if (dimension == 2) {
wf = (1.680351548 * (cos(arg) + 1.)) / hsq;
wfd = -2.639490040 * sin(arg) / (hsq * h);
} else {
wf = 2.051578323 * (cos(arg) + 1.) / (hsq * h);
wfd = -3.222611694 * sin(arg) / (hsq * hsq);
}
}
/*
* compute a normalized smoothing kernel only
*/
static inline void Poly6Kernel(const double hsq, const double h, const double rsq, const int dimension, double &wf) {
double tmp = hsq - rsq;
if (dimension == 2) {
wf = tmp * tmp * tmp / (0.7853981635e0 * hsq * hsq * hsq * hsq);
} else {
wf = tmp * tmp * tmp / (0.6382918409e0 * hsq * hsq * hsq * hsq * h);
}
}
/*
* M4 Prime Kernel
*/
static inline void M4PrimeKernel(const double s, double &wf) {
if (s < 1.0) {
//wf = 1.0 - 2.5 * s * s + (3./2.) * s * s * s;
wf = 1.0 - s * s *(2.5 -1.5 *s);
} else if (s < 2.0) {
//wf = 0.5 * (1.0 - s) * ((2.0 - s) * (2.0 - s));
wf = 2.0 + (-4.0 + (2.5 - 0.5 * s)*s)*s;
} else {
wf = 0.0;
}
}
}
#endif /* SMD_KERNEL_FUNCTIONS_H_ */

View File

@ -0,0 +1,476 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <iostream>
#include "math_special.h"
#include "stdio.h"
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/SVD>
using namespace LAMMPS_NS::MathSpecial;
using namespace std;
using namespace Eigen;
#define MIN(A,B) ((A) < (B) ? (A) : (B))
#define MAX(A,B) ((A) > (B) ? (A) : (B))
/* ----------------------------------------------------------------------
linear EOS for use with linear elasticity
input: initial pressure pInitial, isotropic part of the strain rate d, time-step dt
output: final pressure pFinal, pressure rate p_rate
------------------------------------------------------------------------- */
void LinearEOS(double lambda, double pInitial, double d, double dt, double &pFinal, double &p_rate) {
/*
* pressure rate
*/
p_rate = lambda * d;
pFinal = pInitial + dt * p_rate; // increment pressure using pressure rate
//cout << "hurz" << endl;
}
/* ----------------------------------------------------------------------
shock EOS
input:
current density rho
reference density rho0
current energy density e
reference energy density e0
reference speed of sound c0
shock Hugoniot parameter S
Grueneisen parameter Gamma
initial pressure pInitial
time step dt
output:
pressure rate p_rate
final pressure pFinal
------------------------------------------------------------------------- */
void ShockEOS(double rho, double rho0, double e, double e0, double c0, double S, double Gamma, double pInitial, double dt,
double &pFinal, double &p_rate) {
double mu = rho / rho0 - 1.0;
double pH = rho0 * square(c0) * mu * (1.0 + mu) / square(1.0 - (S - 1.0) * mu);
pFinal = (pH + rho * Gamma * (e - e0));
//printf("shock EOS: rho = %g, rho0 = %g, Gamma=%f, c0=%f, S=%f, e=%f, e0=%f\n", rho, rho0, Gamma, c0, S, e, e0);
//printf("pFinal = %f\n", pFinal);
p_rate = (pFinal - pInitial) / dt;
}
/* ----------------------------------------------------------------------
polynomial EOS
input:
current density rho
reference density rho0
coefficients 0 .. 6
initial pressure pInitial
time step dt
output:
pressure rate p_rate
final pressure pFinal
------------------------------------------------------------------------- */
void polynomialEOS(double rho, double rho0, double e, double C0, double C1, double C2, double C3, double C4, double C5, double C6,
double pInitial, double dt, double &pFinal, double &p_rate) {
double mu = rho / rho0 - 1.0;
if (mu > 0.0) {
pFinal = C0 + C1 * mu + C2 * mu * mu + C3 * mu * mu * mu; // + (C4 + C5 * mu + C6 * mu * mu) * e;
} else {
pFinal = C0 + C1 * mu + C3 * mu * mu * mu; // + (C4 + C5 * mu) * e;
}
pFinal = -pFinal; // we want the mean stress, not the pressure.
//printf("pFinal = %f\n", pFinal);
p_rate = (pFinal - pInitial) / dt;
}
/* ----------------------------------------------------------------------
Tait EOS based on current density vs. reference density.
input: (1) reference sound speed
(2) equilibrium mass density
(3) current mass density
output:(1) pressure
(2) current speed of sound
------------------------------------------------------------------------- */
void TaitEOS_density(const double exponent, const double c0_reference, const double rho_reference, const double rho_current,
double &pressure, double &sound_speed) {
double B = rho_reference * c0_reference * c0_reference / exponent;
double tmp = pow(rho_current / rho_reference, exponent);
pressure = B * (tmp - 1.0);
double bulk_modulus = B * tmp * exponent; // computed as rho * d(pressure)/d(rho)
sound_speed = sqrt(bulk_modulus / rho_current);
// if (fabs(pressure) > 0.01) {
// printf("tmp = %f, press=%f, K=%f\n", tmp, pressure, bulk_modulus);
// }
}
/* ----------------------------------------------------------------------
perfect gas EOS
input: gamma -- adiabatic index (ratio of specific heats)
J -- determinant of deformation gradient
volume0 -- reference configuration volume of particle
energy -- energy of particle
pInitial -- initial pressure of the particle
d -- isotropic part of the strain rate tensor,
dt -- time-step size
output: final pressure pFinal, pressure rate p_rate
------------------------------------------------------------------------- */
void PerfectGasEOS(const double gamma, const double vol, const double mass, const double energy, double &pFinal, double &c0) {
/*
* perfect gas EOS is p = (gamma - 1) rho e
*/
if (energy > 0.0) {
pFinal = (1.0 - gamma) * energy / vol;
//printf("gamma = %f, vol%f, e=%g ==> p=%g\n", gamma, vol, energy, *pFinal__/1.0e-9);
c0 = sqrt((gamma - 1.0) * energy / mass);
} else {
pFinal = c0 = 0.0;
}
}
/* ----------------------------------------------------------------------
linear strength model for use with linear elasticity
input: lambda, mu : Lame parameters
input: sigmaInitial_dev, d_dev: initial stress deviator, deviatoric part of the strain rate tensor
input: dt: time-step
output: sigmaFinal_dev, sigmaFinal_dev_rate__: final stress deviator and its rate.
------------------------------------------------------------------------- */
void LinearStrength(const double mu, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, const double dt,
Matrix3d &sigmaFinal_dev__, Matrix3d &sigma_dev_rate__) {
/*
* deviatoric rate of unrotated stress
*/
sigma_dev_rate__ = 2.0 * mu * d_dev;
/*
* elastic update to the deviatoric stress
*/
sigmaFinal_dev__ = sigmaInitial_dev + dt * sigma_dev_rate__;
}
/* ----------------------------------------------------------------------
linear strength model for use with linear elasticity
input: lambda, mu : Lame parameters
input: F: deformation gradient
output: total stress tensor, deviator + pressure
------------------------------------------------------------------------- */
//void PairTlsph::LinearStrengthDefgrad(double lambda, double mu, Matrix3d F, Matrix3d *T) {
// Matrix3d E, PK2, eye, sigma, S, tau;
//
// eye.setIdentity();
//
// E = 0.5 * (F * F.transpose() - eye); // strain measure E = 0.5 * (B - I) = 0.5 * (F * F^T - I)
// tau = lambda * E.trace() * eye + 2.0 * mu * E; // Kirchhoff stress, work conjugate to above strain
// sigma = tau / F.determinant(); // convert Kirchhoff stress to Cauchy stress
//
////printf("l=%f, mu=%f, sigma xy = %f\n", lambda, mu, sigma(0,1));
//
//// E = 0.5 * (F.transpose() * F - eye); // Green-Lagrange Strain E = 0.5 * (C - I)
//// S = lambda * E.trace() * eye + 2.0 * mu * Deviator(E); // PK2 stress
//// tau = F * S * F.transpose(); // convert PK2 to Kirchhoff stress
//// sigma = tau / F.determinant();
//
// //*T = sigma;
//
// /*
// * neo-hookean model due to Bonet
// */
//// lambda = mu = 100.0;
//// // left Cauchy-Green Tensor, b = F.F^T
// double J = F.determinant();
// double logJ = log(J);
// Matrix3d b;
// b = F * F.transpose();
//
// sigma = (mu / J) * (b - eye) + (lambda / J) * logJ * eye;
// *T = sigma;
//}
/* ----------------------------------------------------------------------
linear strength model for use with linear elasticity
input: lambda, mu : Lame parameters
input: sigmaInitial_dev, d_dev: initial stress deviator, deviatoric part of the strain rate tensor
input: dt: time-step
output: sigmaFinal_dev, sigmaFinal_dev_rate__: final stress deviator and its rate.
------------------------------------------------------------------------- */
void LinearPlasticStrength(const double G, const double yieldStress, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev,
const double dt, Matrix3d &sigmaFinal_dev__, Matrix3d &sigma_dev_rate__, double &plastic_strain_increment) {
Matrix3d sigmaTrial_dev, dev_rate;
double J2;
/*
* deviatoric rate of unrotated stress
*/
dev_rate = 2.0 * G * d_dev;
/*
* perform a trial elastic update to the deviatoric stress
*/
sigmaTrial_dev = sigmaInitial_dev + dt * dev_rate; // increment stress deviator using deviatoric rate
/*
* check yield condition
*/
J2 = sqrt(3. / 2.) * sigmaTrial_dev.norm();
if (J2 < yieldStress) {
/*
* no yielding has occured.
* final deviatoric stress is trial deviatoric stress
*/
sigma_dev_rate__ = dev_rate;
sigmaFinal_dev__ = sigmaTrial_dev;
plastic_strain_increment = 0.0;
//printf("no yield\n");
} else {
//printf("yiedl\n");
/*
* yielding has occured
*/
plastic_strain_increment = (J2 - yieldStress) / (3.0 * G);
/*
* new deviatoric stress:
* obtain by scaling the trial stress deviator
*/
sigmaFinal_dev__ = (yieldStress / J2) * sigmaTrial_dev;
/*
* new deviatoric stress rate
*/
sigma_dev_rate__ = sigmaFinal_dev__ - sigmaInitial_dev;
//printf("yielding has occured.\n");
}
}
/* ----------------------------------------------------------------------
Johnson Cook Material Strength model
input:
G : shear modulus
cp : heat capacity
espec : energy / mass
A : initial yield stress under quasi-static / room temperature conditions
B : proportionality factor for plastic strain dependency
a : exponent for plastic strain dpendency
C : proportionality factor for logarithmic plastic strain rate dependency
epdot0 : dimensionality factor for plastic strain rate dependency
T : current temperature
T0 : reference (room) temperature
Tmelt : melting temperature
input: sigmaInitial_dev, d_dev: initial stress deviator, deviatoric part of the strain rate tensor
input: dt: time-step
output: sigmaFinal_dev, sigmaFinal_dev_rate__: final stress deviator and its rate.
------------------------------------------------------------------------- */
void JohnsonCookStrength(const double G, const double cp, const double espec, const double A, const double B, const double a,
const double C, const double epdot0, const double T0, const double Tmelt, const double M, const double dt, const double ep,
const double epdot, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev__,
Matrix3d &sigma_dev_rate__, double &plastic_strain_increment) {
Matrix3d sigmaTrial_dev, dev_rate;
double J2, yieldStress;
double deltaT = espec / cp;
double TH = deltaT / (Tmelt - T0);
TH = MAX(TH, 0.0);
double epdot_ratio = epdot / epdot0;
epdot_ratio = MAX(epdot_ratio, 1.0);
//printf("current temperature delta is %f, TH=%f\n", deltaT, TH);
yieldStress = (A + B * pow(ep, a)) * (1.0 + C * log(epdot_ratio)); // * (1.0 - pow(TH, M));
/*
* deviatoric rate of unrotated stress
*/
dev_rate = 2.0 * G * d_dev;
/*
* perform a trial elastic update to the deviatoric stress
*/
sigmaTrial_dev = sigmaInitial_dev + dt * dev_rate; // increment stress deviator using deviatoric rate
/*
* check yield condition
*/
J2 = sqrt(3. / 2.) * sigmaTrial_dev.norm();
if (J2 < yieldStress) {
/*
* no yielding has occured.
* final deviatoric stress is trial deviatoric stress
*/
sigma_dev_rate__ = dev_rate;
sigmaFinal_dev__ = sigmaTrial_dev;
plastic_strain_increment = 0.0;
//printf("no yield\n");
} else {
//printf("yiedl\n");
/*
* yielding has occured
*/
plastic_strain_increment = (J2 - yieldStress) / (3.0 * G);
/*
* new deviatoric stress:
* obtain by scaling the trial stress deviator
*/
sigmaFinal_dev__ = (yieldStress / J2) * sigmaTrial_dev;
/*
* new deviatoric stress rate
*/
sigma_dev_rate__ = sigmaFinal_dev__ - sigmaInitial_dev;
//printf("yielding has occured.\n");
}
}
/* ----------------------------------------------------------------------
isotropic maximum strain damage model
input:
current strain
maximum value of allowed principal strain
output:
return value is true if any eigenvalue of the current strain exceeds the allowed principal strain
------------------------------------------------------------------------- */
bool IsotropicMaxStrainDamage(const Matrix3d E, const double maxStrain) {
/*
* compute Eigenvalues of strain matrix
*/
SelfAdjointEigenSolver < Matrix3d > es;
es.compute(E); // compute eigenvalue and eigenvectors of strain
double max_eigenvalue = es.eigenvalues().maxCoeff();
if (max_eigenvalue > maxStrain) {
return true;
} else {
return false;
}
}
/* ----------------------------------------------------------------------
isotropic maximum stress damage model
input:
current stress
maximum value of allowed principal stress
output:
return value is true if any eigenvalue of the current stress exceeds the allowed principal stress
------------------------------------------------------------------------- */
bool IsotropicMaxStressDamage(const Matrix3d S, const double maxStress) {
/*
* compute Eigenvalues of strain matrix
*/
SelfAdjointEigenSolver < Matrix3d > es;
es.compute(S); // compute eigenvalue and eigenvectors of strain
double max_eigenvalue = es.eigenvalues().maxCoeff();
if (max_eigenvalue > maxStress) {
return true;
} else {
return false;
}
}
/* ----------------------------------------------------------------------
Johnson-Cook failure model
input:
output:
------------------------------------------------------------------------- */
double JohnsonCookFailureStrain(const double p, const Matrix3d Sdev, const double d1, const double d2, const double d3,
const double d4, const double epdot0, const double epdot) {
double vm = sqrt(3. / 2.) * Sdev.norm(); // von-Mises equivalent stress
if (vm < 0.0) {
cout << "this is sdev " << endl << Sdev << endl;
printf("vm=%f < 0.0, surely must be an error\n", vm);
exit(1);
}
// determine stress triaxiality
double triax = p / (vm + 0.01 * fabs(p)); // have softening in denominator to avoid divison by zero
if (triax < 0.0) {
triax = 0.0;
} else if (triax > 3.0) {
triax = 3.0;
}
// Johnson-Cook failure strain, dependence on stress triaxiality
double jc_failure_strain = d1 + d2 * exp(d3 * triax);
// include strain rate dependency if parameter d4 is defined and current plastic strain rate exceeds reference strain rate
if (d4 > 0.0) { //
if (epdot > epdot0) {
double epdot_ratio = epdot / epdot0;
jc_failure_strain *= (1.0 + d4 * log(epdot_ratio));
//printf("epsdot=%f, epsdot0=%f, factor = %f\n", epdot, epdot0, (1.0 + d4 * log(epdot_ratio)));
//exit(1);
}
}
return jc_failure_strain;
}

View File

@ -0,0 +1,65 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef SMD_MATERIAL_MODELS_H_
#define SMD_MATERIAL_MODELS_H_
using namespace Eigen;
/*
* EOS models
*/
void LinearEOS(double lambda, double pInitial, double d, double dt, double &pFinal, double &p_rate);
void ShockEOS(double rho, double rho0, double e, double e0, double c0, double S, double Gamma, double pInitial, double dt,
double &pFinal, double &p_rate);
void polynomialEOS(double rho, double rho0, double e, double C0, double C1, double C2, double C3, double C4, double C5, double C6,
double pInitial, double dt, double &pFinal, double &p_rate);
void TaitEOS_density(const double exponent, const double c0_reference, const double rho_reference, const double rho_current,
double &pressure, double &sound_speed);
void PerfectGasEOS(const double gamma, const double vol, const double mass, const double energy, double &pFinal__, double &c0);
/*
* Material strength models
*/
void LinearStrength(const double mu, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, const double dt,
Matrix3d &sigmaFinal_dev__, Matrix3d &sigma_dev_rate__);
void LinearPlasticStrength(const double G, const double yieldStress, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev,
const double dt, Matrix3d &sigmaFinal_dev__, Matrix3d &sigma_dev_rate__, double &plastic_strain_increment);
void JohnsonCookStrength(const double G, const double cp, const double espec, const double A, const double B, const double a,
const double C, const double epdot0, const double T0, const double Tmelt, const double M, const double dt, const double ep,
const double epdot, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev__,
Matrix3d &sigma_dev_rate__, double &plastic_strain_increment);
/*
* Damage models
*/
bool IsotropicMaxStrainDamage(const Matrix3d E, const double maxStrain);
bool IsotropicMaxStressDamage(const Matrix3d E, const double maxStrain);
double JohnsonCookFailureStrain(const double p, const Matrix3d Sdev, const double d1, const double d2, const double d3,
const double d4, const double epdot0, const double epdot);
#endif /* SMD_MATERIAL_MODELS_H_ */

287
src/USER-SMD/smd_math.h Normal file
View File

@ -0,0 +1,287 @@
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
//test
#ifndef SMD_MATH_H_
#define SMD_MATH_H_
#include <Eigen/Eigen>
#include <iostream>
using namespace Eigen;
using namespace std;
namespace SMD_Math {
static inline void LimitDoubleMagnitude(double &x, const double limit) {
/*
* if |x| exceeds limit, set x to limit with the sign of x
*/
if (fabs(x) > limit) { // limit delVdotDelR to a fraction of speed of sound
x = limit * copysign(1, x);
}
}
/*
* deviator of a tensor
*/
static inline Matrix3d Deviator(const Matrix3d M) {
Matrix3d eye;
eye.setIdentity();
eye *= M.trace() / 3.0;
return M - eye;
}
/*
* Polar Decomposition M = R * T
* where R is a rotation and T a pure translation/stretch matrix.
*
* The decomposition is achieved using SVD, i.e. M = U S V^T,
* where U = R V and S is diagonal.
*
*
* For any physically admissible deformation gradient, the determinant of R must equal +1.
* However, scenerios can arise, where the particles interpenetrate and cause inversion, leading to a determinant of R equal to -1.
* In this case, the inversion direction is heuristically identified with the eigenvector of the smallest entry of S, which should work for most cases.
* The sign of this corresponding eigenvalue is flipped, the original matrix M is recomputed using the flipped S, and the rotation and translation matrices are
* obtained again from an SVD. The rotation should proper now, i.e., det(R) = +1.
*/
static inline bool PolDec(Matrix3d M, Matrix3d &R, Matrix3d &T, bool scaleF) {
JacobiSVD<Matrix3d> svd(M, ComputeFullU | ComputeFullV); // SVD(A) = U S V*
Vector3d S_eigenvalues = svd.singularValues();
Matrix3d S = svd.singularValues().asDiagonal();
Matrix3d U = svd.matrixU();
Matrix3d V = svd.matrixV();
Matrix3d eye;
eye.setIdentity();
// now do polar decomposition into M = R * T, where R is rotation
// and T is translation matrix
R = U * V.transpose();
T = V * S * V.transpose();
if (R.determinant() < 0.0) { // this is an improper rotation
// identify the smallest entry in S and flip its sign
int imin;
S_eigenvalues.minCoeff(&imin);
S(imin, imin) *= -1.0;
R = M * V * S.inverse() * V.transpose(); // recompute R using flipped stretch eigenvalues
}
/*
* scale S to avoid small principal strains
*/
if (scaleF) {
double min = 0.3; // 0.3^2 = 0.09, should suffice for most problems
double max = 2.0;
for (int i = 0; i < 3; i++) {
if (S(i, i) < min) {
S(i, i) = min;
} else if (S(i, i) > max) {
S(i, i) = max;
}
}
T = V * S * V.transpose();
}
if (R.determinant() > 0.0) {
return true;
} else {
return false;
}
}
/*
* Pseudo-inverse via SVD
*/
static inline void pseudo_inverse_SVD(Matrix3d &M) {
//JacobiSVD < Matrix3d > svd(M, ComputeFullU | ComputeFullV);
JacobiSVD<Matrix3d> svd(M, ComputeFullU); // one Eigevector base is sufficient because matrix is square and symmetric
Vector3d singularValuesInv;
Vector3d singularValues = svd.singularValues();
//cout << "Here is the matrix V:" << endl << V * singularValues.asDiagonal() * U << endl;
//cout << "Its singular values are:" << endl << singularValues << endl;
double pinvtoler = 1.0e-16; // 2d machining example goes unstable if this value is increased (1.0e-16).
for (int row = 0; row < 3; row++) {
if (singularValues(row) > pinvtoler) {
singularValuesInv(row) = 1.0 / singularValues(row);
} else {
singularValuesInv(row) = 0.0;
}
}
M = svd.matrixU() * singularValuesInv.asDiagonal() * svd.matrixU().transpose();
// JacobiSVD < Matrix3d > svd(M, ComputeFullU | ComputeFullV);
//
// Vector3d singularValuesInv;
// Vector3d singularValues = svd.singularValues();
//
// //cout << "Here is the matrix V:" << endl << V * singularValues.asDiagonal() * U << endl;
// //cout << "Its singular values are:" << endl << singularValues << endl;
//
// double pinvtoler = 1.0e-16; // 2d machining example goes unstable if this value is increased (1.0e-16).
// for (int row = 0; row < 3; row++) {
// if (singularValues(row) > pinvtoler) {
// singularValuesInv(row) = 1.0 / singularValues(row);
// } else {
// singularValuesInv(row) = 0.0;
// }
// }
//
// M = svd.matrixU() * singularValuesInv.asDiagonal() * svd.matrixV().transpose();
}
/*
* test if two matrices are equal
*/
static inline double TestMatricesEqual(Matrix3d A, Matrix3d B, double eps) {
Matrix3d diff;
diff = A - B;
double norm = diff.norm();
if (norm > eps) {
printf("Matrices A and B are not equal! The L2-norm difference is: %g\n", norm);
cout << "Here is matrix A:" << endl << A << endl;
cout << "Here is matrix B:" << endl << B << endl;
}
return norm;
}
/* ----------------------------------------------------------------------
Limit eigenvalues of a matrix to upper and lower bounds.
------------------------------------------------------------------------- */
static inline Matrix3d LimitEigenvalues(Matrix3d S, double limitEigenvalue) {
/*
* compute Eigenvalues of matrix S
*/
SelfAdjointEigenSolver < Matrix3d > es;
es.compute(S);
double max_eigenvalue = es.eigenvalues().maxCoeff();
double min_eigenvalue = es.eigenvalues().minCoeff();
double amax_eigenvalue = fabs(max_eigenvalue);
double amin_eigenvalue = fabs(min_eigenvalue);
if ((amax_eigenvalue > limitEigenvalue) || (amin_eigenvalue > limitEigenvalue)) {
if (amax_eigenvalue > amin_eigenvalue) { // need to scale with max_eigenvalue
double scale = amax_eigenvalue / limitEigenvalue;
Matrix3d V = es.eigenvectors();
Matrix3d S_diag = V.inverse() * S * V; // diagonalized input matrix
S_diag /= scale;
Matrix3d S_scaled = V * S_diag * V.inverse(); // undiagonalize matrix
return S_scaled;
} else { // need to scale using min_eigenvalue
double scale = amin_eigenvalue / limitEigenvalue;
Matrix3d V = es.eigenvectors();
Matrix3d S_diag = V.inverse() * S * V; // diagonalized input matrix
S_diag /= scale;
Matrix3d S_scaled = V * S_diag * V.inverse(); // undiagonalize matrix
return S_scaled;
}
} else { // limiting does not apply
return S;
}
}
static inline bool LimitMinMaxEigenvalues(Matrix3d &S, double min, double max) {
/*
* compute Eigenvalues of matrix S
*/
SelfAdjointEigenSolver < Matrix3d > es;
es.compute(S);
if ((es.eigenvalues().maxCoeff() > max) || (es.eigenvalues().minCoeff() < min)) {
Matrix3d S_diag = es.eigenvalues().asDiagonal();
Matrix3d V = es.eigenvectors();
for (int i = 0; i < 3; i++) {
if (S_diag(i, i) < min) {
//printf("limiting eigenvalue %f --> %f\n", S_diag(i, i), min);
//printf("these are the eigenvalues of U: %f %f %f\n", es.eigenvalues()(0), es.eigenvalues()(1), es.eigenvalues()(2));
S_diag(i, i) = min;
} else if (S_diag(i, i) > max) {
//printf("limiting eigenvalue %f --> %f\n", S_diag(i, i), max);
S_diag(i, i) = max;
}
}
S = V * S_diag * V.inverse(); // undiagonalize matrix
return true;
} else {
return false;
}
}
static inline void reconstruct_rank_deficient_shape_matrix(Matrix3d &K) {
JacobiSVD<Matrix3d> svd(K, ComputeFullU | ComputeFullV);
Vector3d singularValues = svd.singularValues();
for (int i = 0; i < 3; i++) {
if (singularValues(i) < 1.0e-8) {
singularValues(i) = 1.0;
}
}
// int imin;
// double minev = singularValues.minCoeff(&imin);
//
// printf("min eigenvalue=%f has index %d\n", minev, imin);
// Vector3d singularVec = U.col(0).cross(U.col(1));
// cout << "the eigenvalues are " << endl << singularValues << endl;
// cout << "the singular vector is " << endl << singularVec << endl;
//
// // reconstruct original K
//
// singularValues(2) = 1.0;
K = svd.matrixU() * singularValues.asDiagonal() * svd.matrixV().transpose();
//cout << "the reconstructed K is " << endl << K << endl;
//exit(1);
}
/* ----------------------------------------------------------------------
helper functions for crack_exclude
------------------------------------------------------------------------- */
static inline bool IsOnSegment(double xi, double yi, double xj, double yj, double xk, double yk) {
return (xi <= xk || xj <= xk) && (xk <= xi || xk <= xj) && (yi <= yk || yj <= yk) && (yk <= yi || yk <= yj);
}
static inline char ComputeDirection(double xi, double yi, double xj, double yj, double xk, double yk) {
double a = (xk - xi) * (yj - yi);
double b = (xj - xi) * (yk - yi);
return a < b ? -1.0 : a > b ? 1.0 : 0;
}
/** Do line segments (x1, y1)--(x2, y2) and (x3, y3)--(x4, y4) intersect? */
static inline bool DoLineSegmentsIntersect(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4) {
char d1 = ComputeDirection(x3, y3, x4, y4, x1, y1);
char d2 = ComputeDirection(x3, y3, x4, y4, x2, y2);
char d3 = ComputeDirection(x1, y1, x2, y2, x3, y3);
char d4 = ComputeDirection(x1, y1, x2, y2, x4, y4);
return (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) && ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0)))
|| (d1 == 0 && IsOnSegment(x3, y3, x4, y4, x1, y1)) || (d2 == 0 && IsOnSegment(x3, y3, x4, y4, x2, y2))
|| (d3 == 0 && IsOnSegment(x1, y1, x2, y2, x3, y3)) || (d4 == 0 && IsOnSegment(x1, y1, x2, y2, x4, y4));
}
}
#endif /* SMD_MATH_H_ */

View File

@ -43,7 +43,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,ANGMOM,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
MESO_E,MESO_CV,MESO_RHO,INAME,DNAME};
MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,INAME,DNAME};
#define BIG INT_MAX
@ -384,6 +384,24 @@ void Set::command(int narg, char **arg)
set(MESO_RHO);
iarg += 2;
} else if (strcmp(arg[iarg],"smd_mass_density") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->smd_flag)
error->all(FLERR,"Cannot set smd_mass_density for this atom style");
set(SMD_MASS_DENSITY);
iarg += 2;
} else if (strcmp(arg[iarg],"smd_contact_radius") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->smd_flag)
error->all(FLERR,"Cannot set smd_contact_radius for this atom style");
set(SMD_CONTACT_RADIUS);
iarg += 2;
} else if (strstr(arg[iarg],"i_") == arg[iarg]) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@ -571,6 +589,10 @@ void Set::set(int keyword)
else if (keyword == MESO_E) atom->e[i] = dvalue;
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
else if (keyword == SMD_MASS_DENSITY) { // set mass from volume and supplied mass density
atom->rmass[i] = atom->vfrac[i] * dvalue;
}
else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue;
// set shape of ellipsoidal particle