ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

"stringify" and "vectorize" processing of per-atom attributs in AtomVec classes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-14 11:06:10 -04:00
parent 1755d06870
commit b16d48aa41
37 changed files with 1596 additions and 1832 deletions
Download Patch File Download Diff File Expand all files Collapse all files

263
src/atom_vec_tri.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,18 +18,17 @@
#include "error.h"
#include "fix.h"
#include "math_const.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::MY_PI;
#define EPSILON 0.001
static constexpr double EPSILON = 0.001;
/* ---------------------------------------------------------------------- */
@ -58,18 +56,18 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = (char *) "molecule radius rmass omega angmom torque tri";
fields_copy = (char *) "molecule radius rmass omega angmom";
fields_comm = (char *) "";
fields_comm_vel = (char *) "omega angmom";
fields_reverse = (char *) "torque";
fields_border = (char *) "molecule radius rmass";
fields_border_vel = (char *) "molecule radius rmass omega";
fields_exchange = (char *) "molecule radius rmass omega angmom";
fields_restart = (char *) "molecule radius rmass omega angmom";
fields_create = (char *) "molecule radius rmass omega angmom tri";
fields_data_atom = (char *) "id molecule type tri rmass x";
fields_data_vel = (char *) "id v omega angmom";
fields_grow = {"molecule", "radius", "rmass", "omega", "angmom", "torque", "tri"};
fields_copy = {"molecule", "radius", "rmass", "omega", "angmom"};
fields_comm = {};
fields_comm_vel = {"omega", "angmom"};
fields_reverse = {"torque"};
fields_border = {"molecule", "radius", "rmass"};
fields_border_vel = {"molecule", "radius", "rmass", "omega"};
fields_exchange = {"molecule", "radius", "rmass", "omega", "angmom"};
fields_restart = {"molecule", "radius", "rmass", "omega", "angmom"};
fields_create = {"molecule", "radius", "rmass", "omega", "angmom", "tri"};
fields_data_atom = {"id", "molecule", "type", "tri", "rmass", "x"};
fields_data_vel = {"id", "v", "omega", "angmom"};
setup_fields();
}
@ -88,7 +86,7 @@ void AtomVecTri::init()
AtomVec::init();
if (domain->dimension != 3)
error->all(FLERR,"Atom_style tri can only be used in 3d simulations");
error->all(FLERR, "Atom_style tri can only be used in 3d simulations");
}
/* ----------------------------------------------------------------------
@ -112,11 +110,9 @@ void AtomVecTri::grow_pointers()
void AtomVecTri::grow_bonus()
{
nmax_bonus = grow_nmax_bonus(nmax_bonus);
if (nmax_bonus < 0)
error->one(FLERR,"Per-processor system is too big");
if (nmax_bonus < 0) error->one(FLERR, "Per-processor system is too big");
bonus = (Bonus *) memory->srealloc(bonus,nmax_bonus*sizeof(Bonus),
"atom:bonus");
bonus = (Bonus *) memory->srealloc(bonus, nmax_bonus * sizeof(Bonus), "atom:bonus");
}
/* ----------------------------------------------------------------------
@ -129,7 +125,7 @@ void AtomVecTri::copy_bonus(int i, int j, int delflag)
// if deleting atom J via delflag and J has bonus data, then delete it
if (delflag && tri[j] >= 0) {
copy_bonus_all(nlocal_bonus-1,tri[j]);
copy_bonus_all(nlocal_bonus - 1, tri[j]);
nlocal_bonus--;
}
@ -148,7 +144,7 @@ void AtomVecTri::copy_bonus(int i, int j, int delflag)
void AtomVecTri::copy_bonus_all(int i, int j)
{
tri[bonus[i].ilocal] = j;
memcpy(&bonus[j],&bonus[i],sizeof(Bonus));
memcpy(&bonus[j], &bonus[i], sizeof(Bonus));
}
/* ----------------------------------------------------------------------
@ -169,7 +165,7 @@ void AtomVecTri::clear_bonus()
int AtomVecTri::pack_comm_bonus(int n, int *list, double *buf)
{
int i,j,m;
int i, j, m;
double *quat;
m = 0;
@ -191,7 +187,7 @@ int AtomVecTri::pack_comm_bonus(int n, int *list, double *buf)
void AtomVecTri::unpack_comm_bonus(int n, int first, double *buf)
{
int i,m,last;
int i, m, last;
double *quat;
m = 0;
@ -211,13 +207,14 @@ void AtomVecTri::unpack_comm_bonus(int n, int first, double *buf)
int AtomVecTri::pack_border_bonus(int n, int *list, double *buf)
{
int i,j,m;
double *quat,*c1,*c2,*c3,*inertia;
int i, j, m;
double *quat, *c1, *c2, *c3, *inertia;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
if (tri[j] < 0) buf[m++] = ubuf(0).d;
if (tri[j] < 0)
buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
quat = bonus[tri[j]].quat;
@ -251,14 +248,15 @@ int AtomVecTri::pack_border_bonus(int n, int *list, double *buf)
int AtomVecTri::unpack_border_bonus(int n, int first, double *buf)
{
int i,j,m,last;
double *quat,*c1,*c2,*c3,*inertia;
int i, j, m, last;
double *quat, *c1, *c2, *c3, *inertia;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
tri[i] = (int) ubuf(buf[m++]).i;
if (tri[i] == 0) tri[i] = -1;
if (tri[i] == 0)
tri[i] = -1;
else {
j = nlocal_bonus + nghost_bonus;
if (j == nmax_bonus) grow_bonus();
@ -301,7 +299,8 @@ int AtomVecTri::pack_exchange_bonus(int i, double *buf)
{
int m = 0;
if (tri[i] < 0) buf[m++] = ubuf(0).d;
if (tri[i] < 0)
buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
int j = tri[i];
@ -338,7 +337,8 @@ int AtomVecTri::unpack_exchange_bonus(int ilocal, double *buf)
int m = 0;
tri[ilocal] = (int) ubuf(buf[m++]).i;
if (tri[ilocal] == 0) tri[ilocal] = -1;
if (tri[ilocal] == 0)
tri[ilocal] = -1;
else {
if (nlocal_bonus == nmax_bonus) grow_bonus();
double *quat = bonus[nlocal_bonus].quat;
@ -381,8 +381,10 @@ int AtomVecTri::size_restart_bonus()
int n = 0;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (tri[i] >= 0) n += size_restart_bonus_one;
else n++;
if (tri[i] >= 0)
n += size_restart_bonus_one;
else
n++;
}
return n;
@ -396,7 +398,8 @@ int AtomVecTri::pack_restart_bonus(int i, double *buf)
{
int m = 0;
if (tri[i] < 0) buf[m++] = ubuf(0).d;
if (tri[i] < 0)
buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
int j = tri[i];
@ -435,7 +438,8 @@ int AtomVecTri::unpack_restart_bonus(int ilocal, double *buf)
int m = 0;
tri[ilocal] = (int) ubuf(buf[m++]).i;
if (tri[ilocal] == 0) tri[ilocal] = -1;
if (tri[ilocal] == 0)
tri[ilocal] = -1;
else {
if (nlocal_bonus == nmax_bonus) grow_bonus();
double *quat = bonus[nlocal_bonus].quat;
@ -472,52 +476,51 @@ int AtomVecTri::unpack_restart_bonus(int ilocal, double *buf)
void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
{
if (tri[m]) error->one(FLERR,"Assigning tri parameters to non-tri atom");
if (tri[m]) error->one(FLERR, "Assigning tri parameters to non-tri atom");
if (nlocal_bonus == nmax_bonus) grow_bonus();
double c1[3],c2[3],c3[3];
double c1[3], c2[3], c3[3];
int ivalue = 1;
c1[0] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c1[1] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c1[2] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c2[0] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c2[1] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c2[2] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c3[0] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c3[1] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c3[2] = utils::numeric(FLERR,values[ivalue++],true,lmp);
c1[0] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c1[1] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c1[2] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c2[0] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c2[1] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c2[2] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c3[0] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c3[1] = utils::numeric(FLERR, values[ivalue++], true, lmp);
c3[2] = utils::numeric(FLERR, values[ivalue++], true, lmp);
// check for duplicate points
if (c1[0] == c2[0] && c1[1] == c2[1] && c1[2] == c2[2])
error->one(FLERR,"Invalid shape in Triangles section of data file");
error->one(FLERR, "Invalid shape in Triangles section of data file");
if (c1[0] == c3[0] && c1[1] == c3[1] && c1[2] == c3[2])
error->one(FLERR,"Invalid shape in Triangles section of data file");
error->one(FLERR, "Invalid shape in Triangles section of data file");
if (c2[0] == c3[0] && c2[1] == c3[1] && c2[2] == c3[2])
error->one(FLERR,"Invalid shape in Triangles section of data file");
error->one(FLERR, "Invalid shape in Triangles section of data file");
// size = length of one edge
double c2mc1[3],c3mc1[3];
MathExtra::sub3(c2,c1,c2mc1);
MathExtra::sub3(c3,c1,c3mc1);
double size = MAX(MathExtra::len3(c2mc1),MathExtra::len3(c3mc1));
double c2mc1[3], c3mc1[3];
MathExtra::sub3(c2, c1, c2mc1);
MathExtra::sub3(c3, c1, c3mc1);
double size = MAX(MathExtra::len3(c2mc1), MathExtra::len3(c3mc1));
// centroid = 1/3 of sum of vertices
double centroid[3];
centroid[0] = (c1[0]+c2[0]+c3[0]) / 3.0;
centroid[1] = (c1[1]+c2[1]+c3[1]) / 3.0;
centroid[2] = (c1[2]+c2[2]+c3[2]) / 3.0;
centroid[0] = (c1[0] + c2[0] + c3[0]) / 3.0;
centroid[1] = (c1[1] + c2[1] + c3[1]) / 3.0;
centroid[2] = (c1[2] + c2[2] + c3[2]) / 3.0;
double dx = centroid[0] - x[m][0];
double dy = centroid[1] - x[m][1];
double dz = centroid[2] - x[m][2];
double delta = sqrt(dx*dx + dy*dy + dz*dz);
double delta = sqrt(dx * dx + dy * dy + dz * dz);
if (delta/size > EPSILON)
error->one(FLERR,"Inconsistent triangle in data file");
if (delta / size > EPSILON) error->one(FLERR, "Inconsistent triangle in data file");
x[m][0] = centroid[0];
x[m][1] = centroid[1];
@ -528,29 +531,29 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
// tri area = 0.5 len(U x V), where U,V are edge vectors from one vertex
double c4[3];
MathExtra::sub3(c1,centroid,c4);
MathExtra::sub3(c1, centroid, c4);
radius[m] = MathExtra::lensq3(c4);
MathExtra::sub3(c2,centroid,c4);
radius[m] = MAX(radius[m],MathExtra::lensq3(c4));
MathExtra::sub3(c3,centroid,c4);
radius[m] = MAX(radius[m],MathExtra::lensq3(c4));
MathExtra::sub3(c2, centroid, c4);
radius[m] = MAX(radius[m], MathExtra::lensq3(c4));
MathExtra::sub3(c3, centroid, c4);
radius[m] = MAX(radius[m], MathExtra::lensq3(c4));
radius[m] = sqrt(radius[m]);
double norm[3];
MathExtra::cross3(c2mc1,c3mc1,norm);
MathExtra::cross3(c2mc1, c3mc1, norm);
double area = 0.5 * MathExtra::len3(norm);
rmass[m] *= area;
// inertia = inertia tensor of triangle as 6-vector in Voigt ordering
double inertia[6];
MathExtra::inertia_triangle(c1,c2,c3,rmass[m],inertia);
MathExtra::inertia_triangle(c1, c2, c3, rmass[m], inertia);
// diagonalize inertia tensor via Jacobi rotations
// bonus[].inertia = 3 eigenvalues = principal moments of inertia
// evectors and exzy_space = 3 evectors = principal axes of triangle
double tensor[3][3],evectors[3][3];
double tensor[3][3], evectors[3][3];
tensor[0][0] = inertia[0];
tensor[1][1] = inertia[1];
tensor[2][2] = inertia[2];
@ -558,10 +561,10 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
tensor[0][2] = tensor[2][0] = inertia[4];
tensor[0][1] = tensor[1][0] = inertia[5];
int ierror = MathEigen::jacobi3(tensor,bonus[nlocal_bonus].inertia,evectors);
if (ierror) error->one(FLERR,"Insufficient Jacobi rotations for triangle");
int ierror = MathEigen::jacobi3(tensor, bonus[nlocal_bonus].inertia, evectors);
if (ierror) error->one(FLERR, "Insufficient Jacobi rotations for triangle");
double ex_space[3],ey_space[3],ez_space[3];
double ex_space[3], ey_space[3], ez_space[3];
ex_space[0] = evectors[0][0];
ex_space[1] = evectors[1][0];
ex_space[2] = evectors[2][0];
@ -575,26 +578,23 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
// enforce 3 orthogonal vectors as a right-handed coordinate system
// flip 3rd vector if needed
MathExtra::cross3(ex_space,ey_space,norm);
if (MathExtra::dot3(norm,ez_space) < 0.0) MathExtra::negate3(ez_space);
MathExtra::cross3(ex_space, ey_space, norm);
if (MathExtra::dot3(norm, ez_space) < 0.0) MathExtra::negate3(ez_space);
// create initial quaternion
MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus[nlocal_bonus].quat);
MathExtra::exyz_to_q(ex_space, ey_space, ez_space, bonus[nlocal_bonus].quat);
// bonus c1,c2,c3 = displacement of c1,c2,c3 from centroid
// in basis of principal axes
double disp[3];
MathExtra::sub3(c1,centroid,disp);
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
disp,bonus[nlocal_bonus].c1);
MathExtra::sub3(c2,centroid,disp);
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
disp,bonus[nlocal_bonus].c2);
MathExtra::sub3(c3,centroid,disp);
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
disp,bonus[nlocal_bonus].c3);
MathExtra::sub3(c1, centroid, disp);
MathExtra::transpose_matvec(ex_space, ey_space, ez_space, disp, bonus[nlocal_bonus].c1);
MathExtra::sub3(c2, centroid, disp);
MathExtra::transpose_matvec(ex_space, ey_space, ez_space, disp, bonus[nlocal_bonus].c2);
MathExtra::sub3(c3, centroid, disp);
MathExtra::transpose_matvec(ex_space, ey_space, ez_space, disp, bonus[nlocal_bonus].c3);
bonus[nlocal_bonus].ilocal = m;
tri[m] = nlocal_bonus++;
@ -607,7 +607,7 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector<std::string> &values)
double AtomVecTri::memory_usage_bonus()
{
double bytes = 0;
bytes += (double)nmax_bonus*sizeof(Bonus);
bytes += (double) nmax_bonus * sizeof(Bonus);
return bytes;
}
@ -620,7 +620,7 @@ void AtomVecTri::create_atom_post(int ilocal)
{
double radius_one = 0.5;
radius[ilocal] = radius_one;
rmass[ilocal] = 4.0*MY_PI/3.0 * radius_one*radius_one*radius_one;
rmass[ilocal] = 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
tri[ilocal] = -1;
}
@ -632,19 +632,22 @@ void AtomVecTri::create_atom_post(int ilocal)
void AtomVecTri::data_atom_post(int ilocal)
{
tri_flag = tri[ilocal];
if (tri_flag == 0) tri_flag = -1;
else if (tri_flag == 1) tri_flag = 0;
else error->one(FLERR,"Invalid tri flag in Atoms section of data file");
if (tri_flag == 0)
tri_flag = -1;
else if (tri_flag == 1)
tri_flag = 0;
else
error->one(FLERR, "Invalid tri flag in Atoms section of data file");
tri[ilocal] = tri_flag;
if (rmass[ilocal] <= 0.0)
error->one(FLERR,"Invalid density in Atoms section of data file");
if (rmass[ilocal] <= 0.0) error->one(FLERR, "Invalid density in Atoms section of data file");
if (tri_flag < 0) {
double radius_one = 0.5;
radius[ilocal] = radius_one;
rmass[ilocal] *= 4.0*MY_PI/3.0 * radius_one*radius_one*radius_one;
} else radius[ilocal] = 0.0;
rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
} else
radius[ilocal] = 0.0;
omega[ilocal][0] = 0.0;
omega[ilocal][1] = 0.0;
@ -663,17 +666,19 @@ void AtomVecTri::pack_data_pre(int ilocal)
tri_flag = tri[ilocal];
rmass_one = rmass[ilocal];
if (tri_flag < 0) tri[ilocal] = 0;
else tri[ilocal] = 1;
if (tri_flag < 0)
tri[ilocal] = 0;
else
tri[ilocal] = 1;
if (tri_flag < 0) {
double radius_one = radius[ilocal];
rmass[ilocal] /= 4.0*MY_PI/3.0 * radius_one*radius_one*radius_one;
rmass[ilocal] /= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
} else {
double c2mc1[3],c3mc1[3],norm[3];
MathExtra::sub3(bonus[tri_flag].c2,bonus[tri_flag].c1,c2mc1);
MathExtra::sub3(bonus[tri_flag].c3,bonus[tri_flag].c1,c3mc1);
MathExtra::cross3(c2mc1,c3mc1,norm);
double c2mc1[3], c3mc1[3], norm[3];
MathExtra::sub3(bonus[tri_flag].c2, bonus[tri_flag].c1, c2mc1);
MathExtra::sub3(bonus[tri_flag].c3, bonus[tri_flag].c1, c3mc1);
MathExtra::cross3(c2mc1, c3mc1, norm);
double area = 0.5 * MathExtra::len3(norm);
rmass[ilocal] /= area;
}
@ -696,9 +701,9 @@ void AtomVecTri::pack_data_post(int ilocal)
int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
{
int i,j;
double xc,yc,zc;
double dc1[3],dc2[3],dc3[3];
int i, j;
double xc, yc, zc;
double dc1[3], dc2[3], dc3[3];
double p[3][3];
double **x = atom->x;
@ -711,10 +716,10 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
if (buf) {
buf[m++] = ubuf(tag[i]).d;
j = tri[i];
MathExtra::quat_to_mat(bonus[j].quat,p);
MathExtra::matvec(p,bonus[j].c1,dc1);
MathExtra::matvec(p,bonus[j].c2,dc2);
MathExtra::matvec(p,bonus[j].c3,dc3);
MathExtra::quat_to_mat(bonus[j].quat, p);
MathExtra::matvec(p, bonus[j].c1, dc1);
MathExtra::matvec(p, bonus[j].c2, dc2);
MathExtra::matvec(p, bonus[j].c3, dc3);
xc = x[i][0];
yc = x[i][1];
zc = x[i][2];
@ -727,7 +732,8 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
buf[m++] = xc + dc3[0];
buf[m++] = yc + dc3[1];
buf[m++] = zc + dc3[2];
} else m += size_data_bonus;
} else
m += size_data_bonus;
}
return m;
}
@ -740,9 +746,8 @@ void AtomVecTri::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
{
int i = 0;
while (i < n) {
fmt::print(fp,"{} {} {} {} {} {} {} {} {} {}\n", ubuf(buf[i]).i,
buf[i+1],buf[i+2],buf[i+3],buf[i+4],buf[i+5],buf[i+6],
buf[i+7],buf[i+8],buf[i+9]);
fmt::print(fp, "{} {} {} {} {} {} {} {} {} {}\n", ubuf(buf[i]).i, buf[i + 1], buf[i + 2],
buf[i + 3], buf[i + 4], buf[i + 5], buf[i + 6], buf[i + 7], buf[i + 8], buf[i + 9]);
i += size_data_bonus;
}
}
@ -770,24 +775,24 @@ void AtomVecTri::set_equilateral(int i, double size)
quat[1] = 0.0;
quat[2] = 0.0;
quat[3] = 0.0;
c1[0] = -size/2.0;
c1[1] = -sqrt(3.0)/2.0 * size / 3.0;
c1[0] = -size / 2.0;
c1[1] = -sqrt(3.0) / 2.0 * size / 3.0;
c1[2] = 0.0;
c2[0] = size/2.0;
c2[1] = -sqrt(3.0)/2.0 * size / 3.0;
c2[0] = size / 2.0;
c2[1] = -sqrt(3.0) / 2.0 * size / 3.0;
c2[2] = 0.0;
c3[0] = 0.0;
c3[1] = sqrt(3.0)/2.0 * size * 2.0/3.0;
c3[1] = sqrt(3.0) / 2.0 * size * 2.0 / 3.0;
c3[2] = 0.0;
inertia[0] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[1] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[2] = sqrt(3.0)/48.0 * size*size*size*size;
inertia[0] = sqrt(3.0) / 96.0 * size * size * size * size;
inertia[1] = sqrt(3.0) / 96.0 * size * size * size * size;
inertia[2] = sqrt(3.0) / 48.0 * size * size * size * size;
radius[i] = MathExtra::len3(c1);
bonus[nlocal_bonus].ilocal = i;
tri[i] = nlocal_bonus++;
} else if (size == 0.0) {
radius[i] = 0.5;
copy_bonus_all(nlocal_bonus-1,tri[i]);
copy_bonus_all(nlocal_bonus - 1, tri[i]);
nlocal_bonus--;
tri[i] = -1;
} else {
@ -795,18 +800,18 @@ void AtomVecTri::set_equilateral(int i, double size)
double *c2 = bonus[tri[i]].c2;
double *c3 = bonus[tri[i]].c3;
double *inertia = bonus[tri[i]].inertia;
c1[0] = -size/2.0;
c1[1] = -sqrt(3.0)/2.0 * size / 3.0;
c1[0] = -size / 2.0;
c1[1] = -sqrt(3.0) / 2.0 * size / 3.0;
c1[2] = 0.0;
c2[0] = size/2.0;
c2[1] = -sqrt(3.0)/2.0 * size / 3.0;
c2[0] = size / 2.0;
c2[1] = -sqrt(3.0) / 2.0 * size / 3.0;
c2[2] = 0.0;
c3[0] = 0.0;
c3[1] = sqrt(3.0)/2.0 * size * 2.0/3.0;
c3[1] = sqrt(3.0) / 2.0 * size * 2.0 / 3.0;
c3[2] = 0.0;
inertia[0] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[1] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[2] = sqrt(3.0)/48.0 * size*size*size*size;
inertia[0] = sqrt(3.0) / 96.0 * size * size * size * size;
inertia[1] = sqrt(3.0) / 96.0 * size * size * size * size;
inertia[2] = sqrt(3.0) / 48.0 * size * size * size * size;
radius[i] = MathExtra::len3(c1);
}
}